module SS_JPEM_Module

!
! This version of the model implements stochastic choice on a' decisions as a step towards the full
! long term debt model. It also includes coarse / fine grids for decisions / moments and Sunny Day 
! adjustments to control for the extra precautionary savings motive
! 

use kind_mod
use misc_utilities
use omp_lib

implicit none

!--------------------!
! PROGRAM PARAMETERS !
!--------------------!
character(90), parameter :: MasterInDir = 'input/'		! directory with all input files 
character(90), parameter :: MasterOutDir = 'output/'	! directory for all output files 
character(90), parameter :: MasterPanelDir = 'Panel/'	! directory for all output files
character(90) :: InDir, OutDir, PanelDir				! defined in main program based on model version 

integer(ik), parameter :: WReadIn = 1_ik 		! 1 = read in an initial guess of the VALUE FUNCTION from an external file 
integer(ik), parameter :: muReadIn = 1_ik 		! 1 = read in an initial guess of the STATIONARY DISTRIBUTION from an external file 
integer(ik), parameter :: pAReadIn = 1_ik 		! 1 = read in an initial guess of the A' DECISION RULE from an external file 
integer(ik), parameter :: pDReadIn = 1_ik 		! 1 = read in an initial guess of the DEFAULT DECISION RULE from an external file 
integer(ik), parameter :: qReadIn = 1_ik 		! 1 = read in an initial guess of the PRICE FUNCTION from an external file 

integer(ik), parameter :: qMaxIters = 30000_ik 	! max price iterations
integer(ik), parameter :: muMaxIters = 300_ik 	! max stationary distribution iterations 

real(rk), parameter :: realmax = 1.0e+10	! biggest real number (for infinity)
real(rk), parameter :: qTol = 1.0e-3		! tolerance for convergence of price function
real(rk), parameter :: pTol = 1.0e-2 		! tolerance for convergence of decision rules 
real(rk), parameter :: WTol = 1.0e-1 		! tolerance for convergence of value functions 
real(rk), parameter :: muTol = 1.0e-6		! tolerance for convergence of stationary distribution

real(rk), parameter :: WRelax = 1.0_rk 		! weight on new value function for updates
real(rk), parameter :: qRelax = 0.1_rk 		! weight on new pricing function for updates
real(rk), parameter :: pDRelax = qRelax 	! weight on new default decision rule for updates
real(rk), parameter :: pARelax = qRelax 	! weight on new saving decision rule for updates

real(rk), parameter :: gamma_E = 0.5772156649015328606065_rk	! Euler-Mascheroni constant 

!--------------------------------!
! EXTERNALLY ASSIGNED PARAMETERS !
!--------------------------------!
real(rk), parameter :: crra = 3.0_rk  							! crra
real(rk), parameter :: cMin = 0.01_rk							! minimum consumption, LIQUIDATION

real(rk), parameter :: alpha = 0.36_rk							! capital share
real(rk), parameter :: rf_ann = 0.0535_rk						! 5.5% prime rate
real(rk), parameter :: macaulay_0 = 8.32_rk * 4.0_rk / 12.0_rk	! Macaulay duration of debt in months, convert to quarters 
real(rk) :: macaulay = macaulay_0 
real(rk), parameter :: qBar = 1.0_rk 	! NORMALIZATION
real(rk), parameter :: rf = (1.0_rk + rf_ann) ** 0.25_rk - 1.0_rk
real(rk) :: phi 	! decay parameter for long-term bonds (depends on risk-free rate, so determine in eqm loop)
real(rk) :: kappa 	! coupon, used to normalize the risk-free price to 1, also depends on risk-free rate so determine in eqm
real(rk) :: wage 

real(rk), parameter :: psi = 0.0038_rk				! fixed bankruptcy filing cost
real(rk), parameter :: theta_7 = 1.0_rk / 28.0_rk	! probability of regaining credit access

real(rk), parameter :: xi_7_0 = 0.1_rk 				! statutory recovery rate in LIQUIDATION
real(rk) :: xi_7 = xi_7_0

real(rk), parameter :: lambdaBar = 0.4_rk 			! maximal garnishment rate
real(rk), parameter :: xi_13 = 0.2_rk 				! recovery rate rate in PAYMENT PLAN (see Wenli's work)
real(rk), parameter :: theta_13 = 1.0_rk / 20.0_rk	! 1 / avg length of payment plan

integer(ik), parameter :: nXi = 100_ik 	! number of points for discrete approximation to the PDF
real(rk), parameter :: xiMin = 0.0_rk 
real(rk) :: xiMax 
real(rk), allocatable, dimension(:) :: xiGrid, xiProb


!-------!
! GRIDS !
!-------!
! garnishment if under Ch 13 payment plan 
integer(ik), parameter :: nl = 5_ik
real(rk), allocatable, dimension(:) :: lambdaGrid

! QUARTERLY earnings process 
integer(ik), parameter :: ne1 = 5_ik 		! permanent component
integer(ik), parameter :: ne2 = 5_ik		! persistent component 
integer(ik), parameter :: ne3 = 10_ik		! transitory component

real(rk), parameter :: nsd1 = 2.5_rk 		! number of SDs to span for permanent component 
real(rk), parameter :: nsd2 = 2.5_rk 		! number of SDs to span for persistent component 
real(rk), parameter :: nsd3 = 2.5_rk 		! number of SDs to span for transitory component 

real(rk), parameter :: sigma1 = 0.448_rk	! SD, permanent 
real(rk), parameter :: rho2 = 0.984_rk		! persistence, persistent
real(rk), parameter :: sigma2 = 0.080_rk	! SD, persistent
real(rk), parameter :: sigma3 = 0.658_rk	! SD, transitory

real(rk), allocatable, dimension(:) :: e1Grid, e1Share
real(rk), allocatable, dimension(:) :: e2Stat, e2Grid_Orig
real(rk), allocatable, dimension(:, :)  :: e2Trans, cum_e2Trans
real(rk), allocatable, dimension(:,:) :: yBarMat 	! (ne1 x ne2) matrix of expected presented discounted value of labor income over life of Ch 13 bankruptcy
real(rk), allocatable, dimension(:) :: e3Grid_Orig, e3Prob, cum_e3Prob
real(rk) :: AggregateLabor
integer(ik), allocatable, dimension(:) :: e3Map, e2Map 

! exogenous idiosyncratic states combined into a single grid / transition function for simplicity
integer(ik), parameter :: ni = ne2 * ne3
real(rk), allocatable, dimension(:) :: e2Grid, e3Grid, iStat
real(rk), allocatable, dimension(:,:) :: exoShare, iTrans
real(rk), allocatable, dimension(:,:,:) :: lpMat, yMat 

! assets
integer(ik), parameter :: naNeg = 35_ik, naPos = 80_ik				! number of negative and positive asset grid points 
integer(ik), parameter :: na = naNeg + naPos						! total number of asset grid points 
integer(ik), parameter :: aZeroInd = naNeg + 1_ik					! location of a = 0 on the asset grid 
real(rk), parameter :: amin_neg = -1.75_rk, amax_neg = -0.001_rk		! bounds of negative range of the asset grid 
real(rk), parameter :: amin_pos = 0.0_rk, amax_pos = 50.0_rk		! bounds of positive range of the asset grid 
real(rk), allocatable, dimension(:) :: aGridNeg, aGridPos, aGrid	! asset grid, plus negative and positive regions.

!----------------------------------------------!
! INTERNALLY CALIBRATED PARAMETERS AND TARGETS !
!----------------------------------------------!
integer(ik), parameter :: ncal = 10_ik		! number of internally calibrated parameters 

! baseline
real(rk), parameter :: beta_0 = 0.95_rk					! beta 
real(rk), parameter :: stigma_7_0 = 0.5_rk 				! utility cost of LIQUIDATION
real(rk), parameter :: stigma_13_0 = 2.0_rk 			! utility cost of PAYMENT PLAN
real(rk), parameter :: stigma_R_0 = 0.15_rk				! utility cost of DELIQUENCY 
real(rk), parameter :: ev_nu_0 = 0.008_rk	 			! extreme value scale, DEFAULT V NO DEFAULT 
real(rk), parameter :: ev_zeta_0 = 0.005_rk				! extreme value scale, DEFAULT TYPE 
real(rk), parameter :: ev_eta_0 = 0.25_rk				! extreme value scale, BORROWING / SAVING (small, purely numerical)
real(rk), parameter :: IntermediationCost_0 = 0.059_rk	! Intermediation cost per unit
real(rk), parameter :: LoanFixedCost_0 = 0.00000001_rk	! fixed cost of lending (0 in calibrated model)
real(rk), parameter :: alpha_R_0 = 0.0000001_rk			! distribution parameter for exponential discharge in DQ distribution (0 in calibrated model)

real(rk) :: beta = beta_0
real(rk) :: stigma_7 = stigma_7_0
real(rk) :: stigma_13 = stigma_13_0
real(rk) :: stigma_R = stigma_R_0
real(rk) :: ev_nu = ev_nu_0
real(rk) :: ev_zeta = ev_zeta_0
real(rk) :: ev_eta = ev_eta_0
real(rk) :: IntermediationCost = IntermediationCost_0
real(rk) :: LoanFixedCost = LoanFixedCost_0
real(rk) :: alpha_R = alpha_R_0

real(rk), allocatable, dimension(:) :: InternalParameters
! upper and lower bounds of parameter space to search
real(rk), parameter :: beta_LB = 0.9_rk, beta_UB = 0.999_rk
real(rk), parameter :: stigma_7_LB = -7.0_rk, stigma_7_UB = 10.0_rk 
real(rk), parameter :: stigma_13_LB = -10.0_rk, stigma_13_UB = 10.0_rk
real(rk), parameter :: stigma_R_LB = -10.0_rk, stigma_R_UB = 10.0_rk
real(rk), parameter :: ev_nu_LB = 0.000001_rk, ev_nu_UB = 5.0_rk 
real(rk), parameter :: ev_zeta_LB = 0.000001_rk, ev_zeta_UB = 5.0_rk 
real(rk), parameter :: ev_eta_LB = 0.000001_rk, ev_eta_UB = 10.0_rk 
real(rk), parameter :: IntermediationCost_LB = 0.005_rk, IntermediationCost_UB = 0.06_rk
real(rk), parameter :: LoanFixedCost_LB = 0.0_rk, LoanFixedCost_UB = LoanFixedCost_0 + 0.001_rk
real(rk), parameter :: alpha_R_LB = 0.0_rk, alpha_R_UB = 0.4_rk

real(rk), allocatable, dimension(:) :: InternalParameters_LB, InternalParameters_UB
real(rk), allocatable, dimension(:) :: InternalParametersOverReal, ScaledInternalParameters, ExternalParameters

!--------------------!
! ESTIMATION TARGETS !
!--------------------!
integer(ik), parameter :: ntarg = 9_ik					! number of calibration targets 
! empirical targets 
! -- standards 
real(rk), parameter :: Data_BankruptcyRate = 0.404_rk			! bankruptcy rate, pp
real(rk), parameter :: Data_Ch7Share = 57.3_rk					! share of all bankruptcies that are Ch 7 (US Courts, 2023, https://www.uscourts.gov/news/2023/10/26/bankruptcy-filings-rise-13-percent)
real(rk), parameter :: Data_RecoveryRate = 16.0_rk		 		! recovery rate conditional on default
real(rk), parameter :: Data_ChargeoffRate = 3.7_rk		 		! charge-off rate on loans, pp
real(rk), parameter :: Data_FractionInDebt = 11.7_rk			! fraction in debt, pp 
real(rk), parameter :: Data_DebtToIncomeRatio = 4.3_rk			! debt to income ratio, pp 
real(rk), parameter :: Data_AverageLoanSpread = 15.3_rk			! average interest rate spread on loans, annualized pp 
real(rk), parameter :: Data_SpreadIntercept = 14.75_rk 			! intercept of linear fit to pricing schedule (see Figure 2 panel (a) from revision)
real(rk), parameter :: Data_SpreadSlope = 0.12_rk 				! slope of linear fit to pricing schedule 
real(rk), parameter :: Data_KYRatio = 2.2_rk 
! weights on each moment 
real(rk), parameter :: Weight_BankruptcyRate = 1.0_rk
real(rk), parameter :: Weight_Ch7Share = 1.0_rk
real(rk), parameter :: Weight_ChargeoffRate = 1.0_rk
real(rk), parameter :: Weight_RecoveryRate = 1.0_rk
real(rk), parameter :: Weight_FractionInDebt = 1.0_rk 
real(rk), parameter :: Weight_DebtToIncomeRatio = 1.0_rk
real(rk), parameter :: Weight_AverageLoanSpread = 2.0_rk
real(rk), parameter :: Weight_SpreadIntercept = 2.0_rk
real(rk), parameter :: Weight_SpreadSlope = 2.0_rk

real(rk), allocatable, dimension(:)	:: Data_TargetMoments, TargetWeights, Model_TargetMoments, FinalModelMoments
real(rk) :: FractionInDebt, DebtToIncome, DebtToIncome_ap, ChargeoffRate, DefaultRate, Ch7Rate, Ch13Rate, RestructuringRate, &
	& AverageSpread, BankruptcyRate, SuboptimalBankruptcyShare, Ch7Share, SDSpread, RecoveryRate, SpreadIntercept, SpreadSlope


!----------------!
! MODEL VERSIONS !
!----------------!
! version codes:
! 	0: Baseline (estimated parameters)
! 	1: Short Term Debt
! 	2: Low Intermediation Cost 
!   3: Fixed Recovery 1: Ch 7 only 
!   4: Fixed Recovery 2: Ch 7 and Ch 13
integer(ik), parameter :: ModelVersion = 4_ik

real(rk), parameter :: macaulay_ShortTerm = 1.0_rk 
real(rk), parameter :: IntermediationCost_NoIntCost = IntermediationCost_0 * 0.8_rk
real(rk), parameter :: xi_7_FixedRecovery = xi_7_0
real(rk), parameter :: stigma_7_FixedRecovery = stigma_7_0

character(90), parameter :: BaseDir = 'Baseline/'
character(90), parameter :: ShortTermDir = 'ShortTerm/'
character(90), parameter :: LowIntCostDir = 'LowIntCost/'
character(90), parameter :: FixedRecoveryDir = 'FixedRecovery/'
character(90), parameter :: BKOnlyDir = 'BKOnly/'

character(90) :: CaseDir = BaseDir	! folder with inputs selected based on the case of the model being run. 


!--------------------------!
! ENDOGENOUS MODEL OBJECTS !
!--------------------------!
! VALUE FUNCTIONS
real(rk), allocatable, dimension(:,:,:) :: W0, W7		! B-O-P expected value function, no flag and Ch 7
real(rk), allocatable, dimension(:,:,:,:) :: W13		! B-O-P expected value function, Ch 13
real(rk), allocatable, dimension(:,:,:) :: WN0, WD0		! expected value of no default and default, respectively
real(rk), allocatable, dimension(:,:,:) :: J7, J13, JR	! value of each default type

! DECISION RULES and related objects 

real(rk), allocatable, dimension(:,:,:) :: pD, pD7, pD13, pDR	! DEFAULT and DEFAULT TYPE decision rules
real(rk), allocatable, dimension(:,:,:,:) :: pA0, pA7			! BORROW / SAVE DECISION rule, no flag and Ch 7
real(rk), allocatable, dimension(:,:) :: pA0_GoDQ				! transition probability based on xi distribution for going delinquent
real(rk), allocatable, dimension(:,:,:,:,:) :: pA13				! BORROW / SAVE DECISION rule, Ch 13
real(rk), allocatable, dimension(:,:,:,:) :: NashR 				! probability of each possible Nash bargaining outcome in a restructuring

! consumption and sunny day weights associated with all non-default actions 
real(rk), allocatable, dimension(:,:,:,:) :: CA0, CA7, dA0, dA7	! no flag and Ch 7
real(rk), allocatable, dimension(:,:,:,:,:) :: CA13, dA13		! Ch 13

! consumption, recovery, and recovery rates associated with all default actions 
real(rk), allocatable, dimension(:,:,:,:) :: CD7, xiTilde7, xi7, xiTilde13, xi13, lpwMat, lambda13
real(rk), allocatable, dimension(:,:,:) :: CD13, xiTildeR, xiR
integer(ik), allocatable, dimension(:,:,:) :: apiMin0, apiMax0, apiMax7, MustDefault0, CantSave0
integer(ik), allocatable, dimension(:,:,:,:) :: apiMax13, lpiMat


! PRICES
real(rk), allocatable, dimension(:,:,:,:) :: q, spreads, rates		! savings / debt prices (na, na, ni, ne1)
real(rk), allocatable, dimension(:,:,:) :: XD, XN

! DISTRIBUTION 
real(rk), allocatable, dimension(:,:,:) :: mu0, mu7
real(rk), allocatable, dimension(:,:,:,:) :: mu13

! FOR MAPPING PRICES TO data
integer(ik) :: nAhead = 4_ik	! number of periods ahead to compute default probabilities 
real(rk), allocatable, dimension(:,:,:) :: ExpDefProb1, Cond_DefProb1	! "D_1" object from writeup 
real(rk), allocatable, dimension(:,:,:,:) :: ExpDefProb, Cond_DefProb	! "D_n" object from writeup 

integer(ik) :: nPartition = 40_ik	! number of partitions of default probability space 
real(rk) :: pDMin = 0.005_rk, pDMax = 0.2_rk 
real(rk), allocatable, dimension(:) :: pdPartition, pDPartitionPlot	! set of upper bounds of default intervals, and centroids of default intervals for plotting
real(rk), allocatable, dimension(:) :: AvgBinSpread, StDevBinSpread, MassPartition
integer(ik) :: LoanWeight = 1_ik		! set to 1 for balance as opposed to frequency weighting 


!---------------------!
! SIMULATION OF PANEL !
!---------------------!
integer(ik), parameter :: nSimPanel = 100000_ik		! number of individuals to simulate for the panel 
integer(ik), parameter :: tSimPanel = 12_ik			! number of periods to simulate each individual in the panel for 
real(rk), allocatable, dimension(:,:) :: cum_iTrans
real(rk), allocatable, dimension(:) :: cum_xiProb
real(rk), allocatable, dimension(:,:,:,:) :: cum_pA0, cum_pA7
real(rk), allocatable, dimension(:,:,:,:,:) :: cum_pA13

integer(ik) :: tmpTest = 1


contains


!================================================================!
!																 !
! SOLVE THE MODEL AND COMPUTE LOSS FUNCTION FOR GIVEN PARAMETERS !
!																 !
!================================================================!
subroutine SolveModelForGivenParameters(CriterionValue, CurrentMoments, CurrentParametersOverReal)

real(rk), intent(out) :: CriterionValue
real(rk), intent(out), dimension(ntarg) :: CurrentMoments
real(rk), intent(in), dimension(ncal)	:: CurrentParametersOverReal

integer(ik)	:: ai, i, e1i, fi, e2i, e3i, e2ip, e3ip, api, ip, j	! indices

real(rk) :: TKVal	! equilibrium prices 
real(rk) :: aw, apw	! internally calibrated parameters 

!----------------------------------------! 
! INITIALIZATIONS FOR CURRENT PARAMETERS !
!----------------------------------------!
call mapbound(InternalParameters, ncal, CurrentParametersOverReal, InternalParameters_LB, InternalParameters_UB)
beta = InternalParameters(1)
stigma_7 = InternalParameters(2)
stigma_13 = InternalParameters(3)
stigma_R = InternalParameters(4)
ev_nu = InternalParameters(5)
ev_zeta = InternalParameters(6)
ev_eta = InternalParameters(7)
IntermediationCost = InternalParameters(8)
LoanFixedCost = InternalParameters(9)
alpha_R = InternalParameters(10)


!----------------------------------------------------!
!  FOR SINGLE RUNS, IMPOSE THE CURRENT MODEL VERSION !
!----------------------------------------------------!
select case (ModelVersion)
	case(1)	
		macaulay = macaulay_ShortTerm
	case(2)
		IntermediationCost = IntermediationCost_NoIntCost
	case(3)
		xi_7 = xi_7_FixedRecovery
		stigma_7 = stigma_7_FixedRecovery
end select 



188 format(a15, f8.5)
print*, ' '
print 188, 'beta       ', beta
print 188, 'stigma_7   ', stigma_7
print 188, 'stigma_13  ', stigma_13
print 188, 'stigma_R   ', stigma_R
print 188, 'ev_nu      ', ev_nu
print 188, 'ev_zeta    ', ev_zeta
print 188, 'ev_eta     ', ev_eta
print 188, 'Int Cost   ', IntermediationCost
print 188, 'Loan FC    ', LoanFixedCost
print 188, 'alpha_R    ', alpha_R
print*, ' '
print 188, 'xi_7       ', xi_7
print*, ' '


! combine (e2, e3) transitions
do e3i = 1,ne3
do e2i = 1,ne2
	i = (e2i - 1_ik) * ne3 + e3i
	e3Map(i) = e3i
	e2Map(i) = e2i
	
	e2Grid(i) = e2Grid_Orig(e2i)
	e3Grid(i) = e3Grid_Orig(e3i)
	do e3ip = 1,ne3
	do e2ip = 1,ne2	
		ip = (e2ip - 1_ik) * ne3 + e3ip		
		iTrans(i, ip) = e2Trans(e2i, e2ip) * e3Prob(e3ip)
	end do
	end do	
end do
end do


! compute the grid and probability for the xi delinquency shocks
xiMax = expInverseCDF(0.99_rk, alpha_R)
xiMax = min(xiMax, 1.0_rk)
print*, xiMin, xiMax
call linspace(xiMin, xiMax, nXi, xiGrid)
xiProb = expPDF(xiGrid, alpha_R)
xiProb = xiProb / sum(xiProb)

pA0_GoDQ = 0.0_rk
do ai = 1,naNeg
do j = 1,nXi 
	call interp_1d(api, apw, (1.0_rk - xiGrid(j)) * aGrid(ai), na, aGrid)
	pA0_GoDQ(api, ai) = pA0_GoDQ(api, ai) + apw * xiProb(j)
	pA0_GoDQ(api+1, ai) = pA0_GoDQ(api+1, ai) + (1.0_rk - apw) * xiProb(j)
end do 
end do

! distribution
if (muReadIn .eq. 0_ik) then 
	iStat = iTrans(1,:)
	do i = 1,100000
		iStat = matmul(iStat, iTrans)
	end do
	call interp_1d(ai, aw, AggregateLabor * (Data_KYRatio ** (1.0_rk / (1.0_rk - alpha))), na, aGrid)
	mu0 = 0.0_rk
	do e1i = 1,ne1 
	do i = 1,ni
		mu0(ai,i,e1i) = aw * iStat(i) * e1Share(e1i)
		mu0(ai+1,i,e1i) = (1.0_rk - aw) * iStat(i) * e1Share(e1i)
	end do
	end do 
	
	mu7 = 0.0_rk 
	mu13 = 0.0_rk 
end if




!----------------!
! MODEL SOLUTION !
!----------------!
phi = (1.0_rk + rf) / macaulay - rf		! implements Macaulay normalization
kappa = phi + rf 						! implements qBar = 1 normalization 
yMat = wage * lpMat 					! converts labor productivity into labor income
if (minval(yMat) .lt. cmin) then
	print*, 'smallest y Value is too small!'
else
	print*, 'y and cMin consistent'
end if 
print*, minval(yMat)
print*, 'rf', rf
print*, 'wage', wage
print*, 'phi', phi
print*, 'kappa', kappa
print*, 'qBar', qBar
print*, ' '


! SOLVE FOR DECISIONS, PRICES, AND DISTRIBUTION
call Chapter7Details() 				! compute the consumption and recovery stuff for Ch 7 bankruptcy 
call Chapter13Details()				! compute the consumption and recovery stuff for Ch 13 bankruptcy 
call ConsumptionAndWeights_Save()	! compute the consumption and sunny day weights for all savings actions 
print*, ' '
print*, ' '
print*, ' '
print*, ' '
print*, ' '
call SolvePE()						! simultaneous iteration of decision rules and price function 
tmpTest = tmpTest+1
call ComputeAggregateCapital(TKVal, mu0, mu7, mu13)

!--------------------------------------------------!
! MOMENT COMPUTATIONS AND LOSS FUNCTION EVALUATION !
!--------------------------------------------------!
rates = 0.0_rk 
where (q .gt. 0.0_rk)
	rates = kappa * (q ** (-1.0_rk)) - phi 	! compute yield-to-maturity rates
end where 
spreads = rates - rf 						! convert rates to spreads

call OnePeriodAheadDefaultProbability(ExpDefProb1)		! find one-period ahead default probability associated with every possible loan.
if (nAhead .gt. 1_ik) then 
	ExpDefProb(:,:,:,1) = ExpDefProb1
	do i = 2,nAhead 
		call NPeriodAheadDefaultProbability(ExpDefProb(:,:,:,i), ExpDefProb(:,:,:,i-1), ExpDefProb1)
	end do
end if
	
	
	
call ComputeSSMoments()
CurrentMoments(1) = ((1.0_rk + BankruptcyRate) ** 4.0_rk - 1.0_rk) * 100.0_rk	! annualize for comparability
CurrentMoments(2) = Ch7Share * 100.0_rk
CurrentMoments(3) = RecoveryRate * 100.0_rk
CurrentMoments(4) = ChargeoffRate * 100.0_rk
CurrentMoments(5) = FractionInDebt * 100.0_rk
if (abs(100.0_rk * DebtToIncome - Data_DebtToIncomeRatio) .lt. abs(100.0_rk * DebtToIncome_ap - Data_DebtToIncomeRatio)) then 
	CurrentMoments(6) = DebtToIncome * 100.0_rk 
else
	CurrentMoments(6) = DebtToIncome_ap * 100.0_rk 
end if 
CurrentMoments(7) = AverageSpread * 100.0_rk
CurrentMoments(8) = SpreadIntercept
CurrentMoments(9) = SpreadSlope

!-----------!
! REPORTING !
!-----------!
401 format(a40, 2a15)
402 format(a40, 2f15.3)
print*, ' '
print 401, ' moment','model','data'
print 402, ' bankruptcy rate                       ', CurrentMoments(1), Data_TargetMoments(1)
print 402, ' Ch 7 share of bankruptcy              ', CurrentMoments(2), Data_TargetMoments(2)
print 402, ' recovery rate                         ', CurrentMoments(3), Data_TargetMoments(3)
print 402, ' chargeoff rate                        ', CurrentMoments(4), Data_TargetMoments(4)
print 402, ' fraction in debt                      ', CurrentMoments(5), Data_TargetMoments(5)
print 402, ' debt to income ratio                  ', CurrentMoments(6), Data_TargetMoments(6)
print 402, ' average interest rate spread on loans ', CurrentMoments(7), Data_TargetMoments(7)
print 402, ' spread intercept                      ', CurrentMoments(8), Data_TargetMoments(8)
print 402, ' spread slope                          ', CurrentMoments(9), Data_TargetMoments(9)
print*, ' '
CriterionValue = 0.0_rk
do i = 1, ntarg
	CriterionValue = CriterionValue + TargetWeights(i) * (100.0_rk * (CurrentMoments(i) / Data_TargetMoments(i) - 1.0_rk)) ** 2.0_rk	! percentage point deviations
end do
print*, ' criterion value                       ', CriterionValue
print*, ' '
print*, ' '
print*, ' '
print*, '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%'
print*, ' '
print*, ' '
print*, ' '
print*, ' '

! save stuff for Table 6
call savingv(CurrentMoments,'Table6Moments_Model',OutDir)
call savingv(Data_TargetMoments,'Table6Moments_Data',OutDir)
call savingv(InternalParameters(1:8),'Table6Parameters_Internal',OutDir)


ExternalParameters(1) = rf_ann
ExternalParameters(2) = crra
ExternalParameters(3) = phi
ExternalParameters(4) = kappa
ExternalParameters(5) = psi
ExternalParameters(6) = theta_7
ExternalParameters(7) = xi_7
ExternalParameters(8) = cMin
ExternalParameters(9) = theta_13
ExternalParameters(10) = xi_13
ExternalParameters(11) = lambdaBar
ExternalParameters(12) = sigma1
ExternalParameters(13) = rho2
ExternalParameters(14) = sigma2
ExternalParameters(15) = sigma3
call savingv(ExternalParameters, 'Table6Parameters_External',OutDir)

if (isnan(CriterionValue) .eq. .true.) then 
	CriterionValue = huge(CriterionValue) 
end if 

! if doing a single run version of the model, save the output 
print*, ' '
print*, 'saving output...'
print*, ' '
call SaveSSOutput()
print*, 'simulating panel...'
print*, ' '
call PanelSimulation(nSimPanel, tSimPanel, ExpDefProb(:,:,:,nAhead))
print*, ' '
print*, 'completed output.'
print*, ' '

end subroutine SolveModelForGivenParameters




















!========================================================================================!
!                                                                                        !
! SOLVE HH DECISIONS, PRICING FUNCTION, AND STATIONARY DISTRIBUTION FOR GIVEN AGGREGATES !
!                                                                                        !
!========================================================================================!
! SIMULTANEOUS Q AND V ITERATION
subroutine SolvePE()

integer(ik) :: Iter 
real(rk) :: Dist, qDist, XDDist, XNDist, PriceDist
real(rk) :: DefDist, pDDist, pD7Dist, pD13Dist, pDRDist
real(rk) :: ADist, pA0Dist, pA7Dist, pA13Dist, pDist
real(rk) :: WDist, W0Dist, WN0Dist, W7Dist, W13Dist

real(rk), dimension(naNeg,na,ni,ne1)	:: Tq		! update of loan price function 

real(rk), dimension(naNeg,ni,ne1)	 :: TpD, TPD7, TpD13, TpDR, TXD, TXN	! update of default and default type decision rules

real(rk), dimension(na,ni,ne1) :: TW0, TWN0			! update of no flag value function (ex ante and conditional on no default)
real(rk), dimension(naPos,ni,ne1) :: TW7			! update of Ch 7 flag value function 
real(rk), dimension(naPos,ni,ne1,nl) :: TW13		! update of Ch 13 flag value function 

real(rk), dimension(na,na,ni,ne1) :: TpA0			! update of no flag borrow / save decision rule 
real(rk), dimension(naPos,naPos,ni,ne1) :: TpA7		! update of Ch 7 flag borrow / save decision rule 
real(rk), dimension(naPos,naPos,ni,ne1,nl) :: TpA13	! update of Ch 13 flag borrow / save decision rule 

integer(ik) :: ai, i, e1i

! INNER LOOP: loan price schedule and decision rules (always single update method here)
Iter = 0_ik         ! reset repayment function iteration counter
qDist = 1.0_rk       ! initialize p function convergence metric

398 format(20a15)
399 format(i15,19f15.5)
print 398, 'q Iter', 'W0Dist', 'WN0Dist', 'W7Dist', 'W13Dist', 'DefDist', 'ADist', 'pDist', 'qDist','XDDist','XNDist'
do while (Iter .lt. qMaxIters .and. (qDist .gt. qTol .or. pDist .gt. pTol))
	Iter = Iter + 1_ik
	
	!---------------------------------------!
	! SOLVE FOR DECISIONS AND UPDATE PRICES !
	!---------------------------------------!
	! compute consumption for debt actions given the currently assumed debt price schedule 
	call Conditional_OnePeriodAheadDefaultProbability(Cond_DefProb1)		! find one-period ahead default probability associated with every possible loan.
	if (nAhead .gt. 1_ik) then 
		Cond_DefProb(:,:,:,1) = Cond_DefProb1
		do i = 2,nAhead 
			call Conditional_NPeriodAheadDefaultProbability(Cond_DefProb(:,:,:,i), Cond_DefProb(:,:,:,i-1), Cond_DefProb1)
		end do
	end if
	call ConsumptionAndWeights_Borrow(q)
	
	! update the value functions and decision rules 
	call vgSingleUpdate(TW0, TWN0, TW7, TW13, TpA0, TpA7, TpA13, TpD, TpD7, TpD13, TpDR, &	! outputs 
		& W0, W7, W13, WN0, pA0, pA7, pA13, pD, pD7, pD13, pDR, q, apiMin0, apiMax0, MustDefault0, apiMax7, apiMax13, &	! inputs
		& CA0, dA0, CA7, dA7, CA13, dA13, CD7, CD13, xiTilde7, xiTilde13, Cond_DefProb(:,:,:,nAhead))
	
	! compute new pricing function implied by current default decision rules
	call UpdatePriceFunction(Tq, TXD, TXN, q, XD, XN, TpD, TpD7, TpD13, TpDR, xi7, xi13, TpA0)
	
	!----------------------!
	! EVALUATE CONVERGENCE !
	!----------------------!
	! value functions
	W0Dist = maxval(abs(W0 - TW0))
	WN0Dist = maxval(abs(WN0 - TWN0))
	W7Dist = maxval(abs(W7 - TW7))
	W13Dist = maxval(abs(W13 - TW13))
	if (ModelVersion .eq. 3_ik) then 
		WDist = maxval((/ W0Dist, WN0Dist, W7Dist /))
	else
		WDist = maxval((/ W0Dist, WN0Dist, W7Dist, W13Dist /))
	end if 
	
	! default decision rules 
	pDDist = maxval(abs(pD - TpD))
	pD7Dist = maxval(abs(pD7 - TpD7))
	pD13Dist = maxval(abs(pD13 - TpD13))
	pDRDist = maxval(abs(pDR - TpDR))
	
	if (ModelVersion .eq. 3_ik) then 
		DefDist = maxval((/ pDDist,	pD7Dist /))
	elseif (ModelVersion .eq. 4_ik) then 
		DefDist = maxval((/ pDDist,	pD7Dist, pD13Dist /))
	else 
		DefDist = maxval((/ pDDist,	pD7Dist, pD13Dist, pDRDist /))
	end if 
	
	
	! borrow / save decision rules 
	pA0Dist = maxval(abs(pA0 - TpA0))
	pA7Dist = maxval(abs(pA7 - TpA7))
	pA13Dist = maxval(abs(pA13 - TpA13))
	
	if (ModelVersion .eq. 3_ik) then 
		ADist = maxval((/ pA0Dist, pA7Dist /))
	else
		ADist = maxval((/ pA0Dist, pA7Dist, pA13Dist /))
	end if 
	
	pDist = max(DefDist, ADist)
	
	! loan prices 
	qDist = maxval(abs(Tq - q))
	XDDist = maxval(abs(TXD - XD))
	XNDist = maxval(abs(TXN - XN))
	PriceDist = maxval((/ qDist, XDDist, XNDist /))
	
	! final convergence metric
	
	if (modulo(Iter, 50_ik) .eq. 0_ik) then 
		print 399, Iter, W0Dist, WN0Dist, W7Dist, W13Dist, DefDist, ADist, pDist, qDist, XDDist, XNDist, maxval(Tq), minval(Tq)
	end if 
	
	! update prices via relaxation and everything else fully
	q = qRelax * Tq + (1.0_rk - qRelax) * q
	XD = qRelax * TXD + (1.0_rk - qRelax) * XD
	XN = qRelax * TXN + (1.0_rk - qRelax) * XN
	
	W0 = WRelax * TW0 + (1.0_rk - WRelax) * W0
	WN0 = WRelax * TWN0 + (1.0_rk - WRelax) * WN0
	W7 = WRelax * TW7 + (1.0_rk - WRelax) * W7
	W13 = WRelax * TW13 + (1.0_rk - WRelax) * W13
	
	pD = pDRelax * TpD + (1.0_rk - pDRelax) * pD
	pD7 = pDRelax * TpD7 + (1.0_rk - pDRelax) * pD7
	pD13 = pDRelax * TpD13 + (1.0_rk - pDRelax) * pD13
	pDR = pDRelax * TpDR + (1.0_rk - pDRelax) * pDR
	
	pA0 = pARelax * TpA0 + (1.0_rk - pARelax) * pA0
	pA7 = pARelax * TpA7 + (1.0_rk - pARelax) * pA7
	pA13 = pARelax * TpA13 + (1.0_rk - pARelax) * pA13
	
end do  ! end of loop over repayment probability function

print*, ' '
print 399, Iter, W0Dist, WN0Dist, W7Dist, W13Dist, DefDist, ADist, pDist, qDist, XDDist, XNDist

call sdSolve(mu0, mu7, mu13, pA0, pA7, pA13, pD, pD7, pD13, pDR, pA0_GoDQ)	! SOLVE FOR STATIONARY DISTRIBUTION
print*, ' '


call NoEVDecisions(W0, W7, W13, WN0, q, apiMin0, apiMax0, MustDefault0, apiMax7, apiMax13, CA0, CA7, CA13, CD7, CD13)

end subroutine SolvePE






!===================================!
!                                   !
! VALUE FUNCTION AND DECISION RULES !
!                                   !
!===================================!
subroutine vgSingleUpdate(TW0, TWN0, TW7, TW13, TpA0, TpA7, TpA13, TpD, TpD7, TpD13, TpDR, &
	& W0_In, W7_In, W13_In, WN0_In, pA0_In, pA7_In, pA13_In, pD_In, pD7_In, pD13_In, pDR_In, &
	& q_In, apiMin0_In, apiMax0_In, MustDefault0_In, apiMax7_In, apiMax13_In, &
	& CA0_In, dA0_In, CA7_In, dA7_In, CA13_In, dA13_In, CD7_In, CD13_In, xiTilde7_In, xiTilde13_In, delta_1)
	
! value functions 
real(rk), intent(in), dimension(na, ni, ne1) :: W0_In, WN0_In	! current ex-ante value and no-default value (for reorganization value), no flag
real(rk), intent(in), dimension(naPos, ni, ne1) :: W7_In		! current ex-ante Ch 7 value function
real(rk), intent(in), dimension(naPos, ni, ne1, nl) :: W13_In	! current ex-ante Ch 13 value function
real(rk), intent(in), dimension(na, na, ni, ne1) :: pA0_In		! current asset decision rules (needed for restructuring)
real(rk), intent(in), dimension(naPos, naPos, ni, ne1) :: pA7_In		! current asset decision rules (needed for restructuring)
real(rk), intent(in), dimension(naPos, naPos, ni, ne1, nl) :: pA13_In		! current asset decision rules (needed for restructuring)
real(rk), intent(in), dimension(naNeg, ni, ne1) :: pD_In, pD7_In, pD13_In, pDR_In		! current asset decision rules (needed for restructuring)

! loan prices 
real(rk), intent(in), dimension(naNeg, na, ni, ne1) :: q_In 	! current price function

! action- and state-specific consumption, sunny day weights, and asset bounds
integer(ik), intent(in), dimension(na, ni, ne1) :: apiMin0_In, apiMax0_In, MustDefault0_In
integer(ik), intent(in), dimension(naPos, ni, ne1) :: apiMax7_In
integer(ik), intent(in), dimension(naPos, ni, ne1, nl) :: apiMax13_In

real(rk), intent(in), dimension(na, na, ni, ne1) :: CA0_In, dA0_In
real(rk), intent(in), dimension(naPos, naPos, ni, ne1) :: CA7_In, dA7_In
real(rk), intent(in), dimension(naPos, naPos, ni, ne1, nl) :: CA13_In, dA13_In 

real(rk), intent(in), dimension(naNeg, ne1, ne2, ne3) :: CD7_In, xiTilde7_In, xiTilde13_In
real(rk), intent(in), dimension(ne1, ne2, ne3) :: CD13_In

real(rk), intent(in), dimension(naNeg, ni, ne1) :: delta_1

! updated value functions 
real(rk), intent(out), dimension(na, ni, ne1) :: TW0, TWN0 		! updated ex-ante value
real(rk), intent(out), dimension(naPos, ni, ne1) :: TW7 		! updated ex-ante Ch 7 value
real(rk), intent(out), dimension(naPos, ni, ne1, nl) :: TW13 	! updated ex-ante Ch 13 value

! updated default decisions 
real(rk), intent(out), dimension(naNeg, ni, ne1) :: TpD, TpD7, TpD13, TpDR 		! updated default probability and ex-ante value

! updated borrow / save decisions 
real(rk), intent(out), dimension(na, na, ni, ne1) :: TpA0		! updated a' choice probability, no flag / no default
real(rk), intent(out), dimension(naPos, naPos, ni, ne1) :: TpA7	! updated a' choice probability, Ch 7
real(rk), intent(out), dimension(naPos, naPos, ni, ne1, nl) :: TpA13	! updated a' choice probability, Ch 7

! TEMPORARIES
integer(ik)	:: ai, e1i, i, apiMin, apiMax, api, aiPos, li, lpi, e2i, e3i, j
real(rk) :: yVal, cDVal, WDVal, WNVal, WMaxVal, zSum, aVal, lpw, J7Val, J13Val, JRVal, JMaxVal, zJSum, WRBarVal, apw
real(rk), dimension(3) :: JVec, zJVec
real(rk), dimension(na,ni,ne1) :: W0Bar  			! expected value next period, no flag today and no default
real(rk), dimension(ni,ne1) :: W0Bar_FileCh7 		! expected value next period, no flag today and Ch 7
real(rk), dimension(naNeg,ni,ne1) :: W0Bar_FileCh13, W0Bar_GoDQ ! expected value next period, no flag today, file Ch 13 and go Delinquent

real(rk), dimension(naPos,ni,ne1) :: W7Bar  		! expected value next period, Ch7 today
real(rk), dimension(naPos,ni,ne1,nl) :: W13Bar  	! expected value next period, Ch13 today

!---------------------------------------------------!
! COMPUTE CONDITIONAL EXPECTATION OF VALUE FUNCTION !
!---------------------------------------------------!
do e1i = 1, ne1
do i = 1, ni
	! no default flag, continuation value associated with filing for Ch7 default today
	W0Bar_FileCh7(i,e1i) = dot_product(iTrans(i,:), W7_In(1,:,e1i))		! go to Ch 7 for sure, integrate over idiosyncratic shocks 
	e2i = e2Map(i)
	e3i = e3Map(i)
	
	!$omp parallel do private(aiPos, lpi, lpw)
	do ai = 1, na
		! no default flag, no default decision 
		W0Bar(ai,i,e1i) = dot_product(iTrans(i,:), W0_In(ai,:,e1i))
		
		if (ai .lt. aZeroInd) then 
			! no default flag, file for Ch 13 today: go to a' = 0, get Ch 13 flag, and get lambda implied by Ch 13 recovery determination procedure
			lpi = lpiMat(ai, e1i, e2i, e3i)
			lpw = lpwMat(ai, e1i, e2i, e3i)
			W0Bar_FileCh13(ai,i,e1i) = dot_product(iTrans(i,:), lpw * W13_In(1,:,e1i,lpi) + (1.0_rk - lpw) * W13_In(1,:,e1i,lpi+1))
		else
			! default flag 
			aiPos = aiPositive(ai)
			W7Bar(aiPos,i,e1i) = dot_product(iTrans(i,:), theta_7 * W0_In(ai,:,e1i) + (1.0_rk - theta_7) * W7_In(aiPos,:,e1i))
			do li = 1,nl
				W13Bar(aiPos,i,e1i,li) = dot_product(iTrans(i,:), theta_13 * W0_In(ai,:,e1i) + (1.0_rk - theta_13) * W13_In(aiPos,:,e1i,li))
			end do 
		end if 
	end do  ! ai, expected value
	!$omp end parallel do
end do  ! i, expected value
end do  ! e1i, expected value

! Compute the conditional expectation associated with going delinquent
do e1i = 1,ne1 
do i = 1,ni 
!$omp parallel do private(WRBarVal, api, apw)
do ai = 1,naNeg
	WRBarVal = 0.0_rk 
	do j = 1,nXi 
		call interp_1d(api, apw, (1.0_rk - xiGrid(j)) * aGrid(ai), na, aGrid)
		WRBarVal = WRBarVal + xiProb(j) * (apw * W0Bar(api, i, e1i) + (1.0_rk - apw) * W0Bar(api+1, i, e1i))
	end do 
	W0Bar_GoDQ(ai,i,e1i) = WRBarVal
end do
!$omp end parallel do
end do
end do 

!-------------------------------!
! APPLY THE FUNCTIONAL OPERATOR !
!-------------------------------!
TpA0 = pA0_In
TpA7 = pA7_In
TpA13 = pA13_In

TpD = pD_In
TpD7 = pD7_In
TpD13 = pD13_In
TpDR = pDR_In

TW0 = W0_In
TWN0 = WN0_In 
TW7 = W7_In
TW13 = W13_In

NashR = 0.0_rk 

do e1i = 1, ne1
do i = 1, ni
	e2i = e2Map(i)
	e3i = e3Map(i) 
	yVal = yMat(e1i,e2i,e3i)
	
	!$omp parallel do private(aVal, WNVal, J7Val, J13Val, JRVal, JVec, JMaxVal, zJVec, zJSum, WDVal, WMaxVal, zSum, aiPos)
	do ai = 1,na
		aVal = aGrid(ai)
		! REPAYMENT, NO FLAG
		if (MustDefault0_In(ai,i,e1i) .eq. 0_ik) then 		! repayment is feasible 
			call SolveBorrowSave(WNVal, TpA0(:,ai,i,e1i), apiMin0_In(ai,i,e1i), apiMax0_In(ai,i,e1i), &
				& CA0_In(:,ai,i,e1i), dA0_In(:,ai,i,e1i), W0Bar(:,i,e1i), na, delta_1(:,i,e1i), ai)
		else 
			WNVal = -realmax
		end if 
		
		! NEGATIVE A RANGE: SOLVE FOR DEFAULT, ONLY FOR NO FLAG
		if (ai .le. naNeg) then 
			! type-specific values (for default actions that are always feasible)
			J13Val = utility(cD13_In(e1i, e2i, e3i)) - stigma_13 + beta * W0Bar_FileCh13(ai, i, e1i)	! Ch 13 
			JRVal = utility(yVal) - stigma_R + beta * W0Bar_GoDQ(ai, i, e1i)							! Delinquency / restructuring
			if (yVal .lt. 1.0_rk / wage) then 																	! adjust for feasibility of Ch 7 given income restrictions
				J7Val = utility(cD7_In(ai, e1i, e2i, e3i)) - stigma_7 + beta * W0Bar_FileCh7(i,e1i)		! Ch 7 
			else 
				J7Val = -realmax
			end if 
			
			if (ModelVersion .ge. 3_ik) then
				JRVal = -realmax
			end if 
			
			if (ModelVersion .eq. 3_ik) then
				J13Val = -realmax
			end if 
			
			JVec = (/ J7Val, J13Val, JRVal /)
			JMaxVal = maxval(JVec)
			zJVec = exp((JVec - JMaxVal) / ev_zeta)
			zJSum = sum(zJVec)
			WDVal = JMaxVal + ev_zeta * log(zJSum) ! - ev_zeta * log(3.0_rk)
			TpD7(ai,i,e1i) = zJVec(1) / zJSum
			TpD13(ai,i,e1i) = zJVec(2) / zJSum
			TpDR(ai,i,e1i) = (1.0_rk - TpD7(ai,i,e1i) - TpD13(ai,i,e1i))

			! SOLVE DEFAULT / NO DEFAULT DECISION
			WMaxVal = max(WDVal, WNVal)
			zSum = exp((WDVal - WMaxVal) / ev_nu) + exp((WNVal - WMaxVal) / ev_nu)
			TW0(ai,i,e1i) = WMaxVal + ev_nu * log(zSum) - ev_nu * log(2.0_rk)
			TpD(ai,i,e1i) = exp((WDVal - WMaxVal) / ev_nu) / zSum
			
			TWN0(ai,i,e1i) = WNVal
			J7(ai,i,e1i) = J7Val 
			J13(ai,i,e1i) = J13Val
			JR(ai,i,e1i) = JRVal
			WD0(ai,i,e1i) = WDVal
		else 
			! IN POSITIVE RANGE, DON'T HAVE TO SOLVE FOR DEFAULT, BUT DO HAVE TO SOLVE THE PROBLEMS FOR OTHER FLAG STATUSES !
			aiPos = aiPositive(ai)
			
			! SOLVE FOR CH 7 AND CH 13 VALUES
			call SolveBorrowSave(TW7(aiPos,i,e1i), TpA7(:,aiPos,i,e1i), 1_ik, apiMax7_In(aiPos,i,e1i), CA7_In(:,aiPos,i,e1i), &
				& dA7_In(:,aiPos,i,e1i), W7Bar(:,i,e1i), naPos, delta_1(:,i,e1i), ai) 
			do li = 1,nl
				call SolveBorrowSave(TW13(aiPos,i,e1i,li), TpA13(:,aiPos,i,e1i,li), 1_ik, apiMax13_In(aiPos,i,e1i,li), CA13_In(:,aiPos,i,e1i,li), &
					& dA13_In(:,aiPos,i,e1i,li), W13Bar(:,i,e1i,li), naPos, delta_1(:,i,e1i), ai)
			end do 
			
			! NO DEFAULT FLAG: DEFAULT IS INFEASIBLE, SO ASSIGN REPAYMENT VALUE
			TW0(ai,i,e1i) = WNVal
			TWN0(ai,i,e1i) = WNVal
		end if 
	end do
	!$omp end parallel do
end do
end do

end subroutine vgSingleUpdate




!----------------------------------------------------------------!
! SOLVE THE BORROWING / SAVING PROBLEM CONDITIONAL ON NO DEFAULT !
!----------------------------------------------------------------!
! given consumption, sunny day weights, continuation value, and range of feasible actions 
! NB: this is the same across default flag statuses, conditional on providing the correct inputs.
subroutine SolveBorrowSave(WVal, ProbVec, apiMin, apiMax, cVec, dVec, WBarVec, nIn, delta_1, ai)
integer(ik), intent(in) :: apiMin, apiMax, nIn, ai
real(rk), intent(in), dimension(nIn) :: cVec, dVec, WBarVec	! consumption values, sunny day weights, and continuation values 
real(rk), intent(in), dimension(naNeg) :: delta_1

real(rk), intent(out) :: WVal	! ex ante value of feasible borrow / save choices given current state
real(rk), intent(out), dimension(nIn) :: ProbVec

real(rk) :: zSum, vMaxVal
real(rk), dimension(nIn) :: zVec, vVec
integer(ik) :: api 
ProbVec = 0.0_rk 
vVec(apiMin:apiMax) = utility(cVec(apiMin:apiMax)) + beta * WBarVec(apiMin:apiMax)			! value of feasible actions 
vMaxVal = maxval(vVec(apiMin:apiMax))														! max value
zVec(apiMin:apiMax) = dVec(apiMin:apiMax) * exp((vVec(apiMin:apiMax) - vMaxVal) / ev_eta)	! vector of exponential terms 
zSum = sum(zVec(apiMin:apiMax))																! sum term of exponential terms
if (zSum .gt. 0.0_rk) then 																	! control for actions being so bad, most germane to indebted no flag case 
	WVal= vMaxVal + ev_eta * log(zSum)														! ex ante value
	ProbVec(apiMin:apiMax) = zVec(apiMin:apiMax) / zSum										! action probabilities conditional on repayment
	ProbVec = ProbVec / sum(ProbVec)														! clean up 
else
	WVal = -realmax 
end if
end subroutine SolveBorrowSave 















!==================================================!
!                                                  !
! DECISION RULES FOR THE DSSY "EPSILON-ZERO" MODEL !
!                                                  !
!==================================================!
subroutine NoEVDecisions(W0_In, W7_In, W13_In, WN0_In, q_In, apiMin0_In, apiMax0_In, MustDefault0_In, apiMax7_In, apiMax13_In, CA0_In, CA7_In, CA13_In, CD7_In, CD13_In)
	
! value functions 
real(rk), intent(in), dimension(na, ni, ne1) :: W0_In, WN0_In	! current ex-ante value and no-default value (for reorganization value), no flag
real(rk), intent(in), dimension(naPos, ni, ne1) :: W7_In		! current ex-ante Ch 7 value function
real(rk), intent(in), dimension(naPos, ni, ne1, nl) :: W13_In	! current ex-ante Ch 13 value function
! loan prices 
real(rk), intent(in), dimension(naNeg, na, ni, ne1) :: q_In 	! current price function

! action- and state-specific consumption, sunny day weights, and asset bounds
integer(ik), intent(in), dimension(na, ni, ne1) :: apiMin0_In, apiMax0_In, MustDefault0_In
integer(ik), intent(in), dimension(naPos, ni, ne1) :: apiMax7_In
integer(ik), intent(in), dimension(naPos, ni, ne1, nl) :: apiMax13_In

real(rk), intent(in), dimension(na, na, ni, ne1) :: CA0_In
real(rk), intent(in), dimension(naPos, naPos, ni, ne1) :: CA7_In
real(rk), intent(in), dimension(naPos, naPos, ni, ne1, nl) :: CA13_In

real(rk), intent(in), dimension(naNeg, ne1, ne2, ne3) :: CD7_In
real(rk), intent(in), dimension(ne1, ne2, ne3) :: CD13_In

! updated decisions 
integer(ik), dimension(naNeg, ni, ne1) :: gD, gD7, gD13, gDR 	! deterministic default and default type decisions
integer(ik), dimension(na, ni, ne1) :: gai0						! deterministic a' choice, no flag / no default
integer(ik), dimension(naPos, ni, ne1) :: gai7					! deterministic a' choice, Ch 7
integer(ik), dimension(naPos, ni, ne1, nl) :: gai13				! deterministic a' choice, Ch 13

! TEMPORARIES
integer(ik)	:: ai, e1i, i, apiMin, apiMax, api, aiPos, li, lpi, e2i, e3i, j
real(rk) :: yVal, cDVal, WDVal, WNVal, J7Val, J13Val, JRVal, lpw, apw, WRBarVal 
real(rk), dimension(na) :: vVec 
real(rk), dimension(na,ni,ne1) :: W0Bar  			! expected value next period, no flag today and no default
real(rk), dimension(ni,ne1) :: W0Bar_FileCh7 		! expected value next period, no flag today and Ch 7
real(rk), dimension(naNeg,ni,ne1) :: W0Bar_FileCh13, W0Bar_GoDQ ! expected value next period, no flag today, file Ch 13 and go Delinquent

real(rk), dimension(naPos,ni,ne1) :: W7Bar  		! expected value next period, Ch7 today
real(rk), dimension(naPos,ni,ne1,nl) :: W13Bar  	! expected value next period, Ch13 today

character(90) :: Dir

!--------------------------------------------------------------------------------!
! COMPUTE CONDITIONAL EXPECTATION OF VALUE FUNCTION -- SAME AS IN BASELINE MODEL !
!--------------------------------------------------------------------------------!
do e1i = 1, ne1
do i = 1, ni
	! no default flag, continuation value associated with filing for Ch7 default today
	W0Bar_FileCh7(i,e1i) = dot_product(iTrans(i,:), W7_In(1,:,e1i))		! go to Ch 7 for sure, integrate over idiosyncratic shocks 
	e2i = e2Map(i)
	e3i = e3Map(i)
	
	!$omp parallel do private(aiPos, lpi, lpw)
	do ai = 1, na
		! no default flag, no default decision 
		W0Bar(ai,i,e1i) = dot_product(iTrans(i,:), W0_In(ai,:,e1i))
		
		if (ai .lt. aZeroInd) then 
			! no default flag, file for Ch 13 today: go to a' = 0, get Ch 13 flag, and get lambda implied by Ch 13 recovery determination procedure
			lpi = lpiMat(ai, e1i, e2i, e3i)
			lpw = lpwMat(ai, e1i, e2i, e3i)
			W0Bar_FileCh13(ai,i,e1i) = dot_product(iTrans(i,:), lpw * W13_In(1,:,e1i,lpi) + (1.0_rk - lpw) * W13_In(1,:,e1i,lpi+1))
		else
			! default flag 
			aiPos = aiPositive(ai)
			W7Bar(aiPos,i,e1i) = dot_product(iTrans(i,:), theta_7 * W0_In(ai,:,e1i) + (1.0_rk - theta_7) * W7_In(aiPos,:,e1i))
			do li = 1,nl
				W13Bar(aiPos,i,e1i,li) = dot_product(iTrans(i,:), theta_13 * W0_In(ai,:,e1i) + (1.0_rk - theta_13) * W13_In(aiPos,:,e1i,li))
			end do 
		end if 
	end do  ! ai, expected value
	!$omp end parallel do
end do  ! i, expected value
end do  ! e1i, expected value

! Compute the conditional expectation associated with going delinquent
do e1i = 1,ne1 
do i = 1,ni 
!$omp parallel do private(WRBarVal, api, apw)
do ai = 1,naNeg
	WRBarVal = 0.0_rk 
	do j = 1,nXi 
		call interp_1d(api, apw, (1.0_rk - xiGrid(j)) * aGrid(ai), na, aGrid)
		WRBarVal = WRBarVal + xiProb(j) * (apw * W0Bar(api, i, e1i) + (1.0_rk - apw) * W0Bar(api+1, i, e1i))
	end do 
	W0Bar_GoDQ(ai,i,e1i) = WRBarVal
end do
!$omp end parallel do
end do
end do 

!--------------------------------------------------------------------------------------------!
! APPLY THE FUNCTIONAL OPERATOR -- MODIFIED FROM BASELINE MODEL TO KILL EXTREME VALUE SHOCKS !
!--------------------------------------------------------------------------------------------!
NashR = 0.0_rk 
gai0 = aZeroInd
gd7 = 0_ik 
gD13 = 0_ik 
gDR = 0_ik 
gD = 0_ik

do e1i = 1, ne1
do i = 1, ni
	e2i = e2Map(i)
	e3i = e3Map(i) 
	yVal = yMat(e1i,e2i,e3i)
	
	!$omp parallel do private(apiMin, apiMax, vVec, WNVal, J13Val, JRVal, J7Val, WDVal, aiPos)
	do ai = 1,na
		! REPAYMENT, NO FLAG
		if (MustDefault0_In(ai,i,e1i) .eq. 0_ik) then 		! repayment is feasible 
			apiMin = apiMin0_In(ai,i,e1i)
			apiMax = apiMax0_in(ai,i,e1i)
			vVec = -realmax
			vVec(apiMin:apiMax) = utility(cA0_In(apiMin:apiMax,ai,i,e1i)) + beta * W0Bar(1:apiMax,i,e1i)
			gai0(ai,i,e1i) = maxloc(vVec(1:apiMax), 1)
			WNVal = maxval(vVec(apiMin:apiMax))
		else 
			WNVal = -realmax
		end if 
		
		! NEGATIVE A RANGE: SOLVE FOR DEFAULT, ONLY FOR NO FLAG
		if (ai .le. naNeg) then 
			! -- type-specific values (for default actions that are always feasible)
			J13Val = utility(cD13_In(e1i, e2i, e3i)) - stigma_13 + beta * W0Bar_FileCh13(ai, i, e1i)
			JRVal = utility(yVal) - stigma_R + beta * W0Bar_GoDQ(ai, i, e1i)

			! -- ex ante default value
			if (yVal / wage .lt. 1.0_rk) then 				! adjust for feasibility of Ch 7 given income restrictions
				J7Val = utility(cD7_In(ai, e1i, e2i, e3i)) - stigma_7 + beta * W0Bar_FileCh7(i,e1i)
				
				if (J7Val .ge. JRVal .and. J7Val .ge. J13Val) then
					gD7(ai, i, e1i) = 1_ik
					WDVal = J7Val
				elseif (J13Val .ge. JRVal) then
					gD13(ai, i, e1i) = 1_ik
					WDVal = J13Val
				else
					gDR(ai, i, e1i) = 1_ik
					WDVal = JRVal 
				end if 
			else 
				if (J13Val .ge. JRVal) then
					gD13(ai, i, e1i) = 1_ik 
					WDVal = J13Val
				else 
					gDR(ai, i, e1i) = 1_ik 
					WDVal = JRVal 
				end if 
			end if 
			
			! SOLVE DEFAULT / NO DEFAULT DECISION
			if (WDVal .gt. WNVal) then 
				gD(ai, i, e1i) = 1_ik
			end if
		else 
			! IN POSITIVE RANGE, DON'T HAVE TO SOLVE FOR DEFAULT, BUT DO HAVE TO SOLVE THE PROBLEMS FOR OTHER FLAG STATUSES
			aiPos = aiPositive(ai)
			
			! Ch 7 flag 
			apiMax = apiMax7_In(aiPos, i, e1i)
			vVec(1:apiMax) = utility(cA7_In(1:apiMax, aiPos, i, e1i)) + beta * W7Bar(1:apiMax, i, e1i)
			gai7(aiPos, i, e1i) = maxloc(vVec(1:apiMax), 1)
						
			! Ch 13 flag
			do li = 1,nl
				apiMax = apiMax13_In(aiPos, i, e1i, li)
				vVec(1:apiMax) = utility(cA13_In(1:apiMax, aiPos, i, e1i, li)) + beta * W13Bar(1:apiMax, i, e1i, li)
				gai13(aiPos, i, e1i, li) = maxloc(vVec(1:apiMax), 1)
			end do 
		end if 
	end do
	!$omp end parallel do
end do
end do



! SAVE 
Dir = trim(OutDir) // trim('NoEV/')
call saving3d(dble(gD),'gD',Dir)
call saving3d(dble(gD7),'gD7',Dir)
call saving3d(dble(gD13),'gD13',Dir)
call saving3d(dble(gDR),'gDR',Dir)

call saving3d(dble(gai0),'gai0',Dir)
call saving3d(dble(gai7),'gai7',Dir)
call saving4d(dble(gai13),'gai13',Dir)

print*, 'gD'
print*, minval(gD), maxval(gD)

end subroutine NoEVDecisions












!==============================================================!
!                                                              !
! COMPUTE PRICING FUNCTION GIVEN DECISION RULES, DISCOUNT RATE !
!                                                              !
!==============================================================!
subroutine UpdatePriceFunction(Tq, TXD, TXN, q_In, XD_In, XN_In, pD_In, pD7_In, pD13_In, pDR_In, xi7_In, xi13_In, pA0_In)
real(rk), intent(in), dimension(naNeg, ni, ne1) :: pD_In, pD7_In, pD13_In, pDR_In	! current default and default type decision rules
real(rk), intent(in), dimension(naNeg, ne1, ne2, ne3) :: xi7_In, xi13_In					! recovery rates by default type and individual state at the time of default
real(rk), intent(in), dimension(na, na, ni, ne1) :: pA0_In		! current borrow / save decision rule
real(rk), intent(in), dimension(naNeg, na, ni, ne1) :: q_In		! current loan price
real(rk), intent(out), dimension(naNeg, na, ni, ne1) :: Tq		! updated loan price 
real(rk), intent(in), dimension(naNeg, ni, ne1) :: XD_In, XN_In	! current recovery value in default
real(rk), intent(out), dimension(naNeg, ni, ne1) :: TXD, TXN	! updated recovery value in default

real(rk) :: TqVal, LoanSize, DVal, NVal, pDVal, xiRVal, apw
integer(ik) :: e1i, i, ai, api, ip, e2i, e3i, j
real(rk), dimension(3) :: xijVec, pDjVec	! default-type-specific recovery and choice probability 
real(rk), dimension(na) :: qVec 

! update XN
do e1i = 1,ne1 
do i = 1,ni
	e2i = e2Map(i)
	e3i = e3Map(i)
!$omp parallel do private(qVec)
do ai = 1,naNeg 
	qVec(1:naNeg) = q_In(:,ai,i,e1i)	! current debt prices 
	qVec(aZeroInd:na) = qBar			! everything else gets the (normalized) savings price 
	TXN(ai, i, e1i) = kappa + (1.0_rk - phi) * dot_product(pA0_In(:,ai,i,e1i), qVec)
end do
!$omp end parallel do 
end do
end do 

! update XD
do e1i = 1,ne1 
do i = 1,ni
	e2i = e2Map(i)
	e3i = e3Map(i)
!$omp parallel do private(pDjVec, xiRVal, api, apw, DVal, NVal, xijVec)
do ai = 1,naNeg 
	! default value 
	pDjVec = (/ pD7_In(ai, i, e1i), pD13_In(ai, i, e1i), pDR_In(ai, i, e1i) /)	! vector of probabilities of default type 
	xiRVal = 0.0_rk 
	do j = 1,nXi
		call interp_1d(api, apw, (1.0_rk - xiGrid(j)) * aGrid(ai), na, aGrid) 
		do ip = 1,ni
			DVal = apw * pD_In(api, ip, e1i) * XD_In(api,ip,e1i) + (1.0_rk - apw) * pD_In(api+1, ip, e1i) * XD_In(api+1,ip,e1i)
			NVal = apw * (1.0_rk - pD_In(api, ip, e1i)) * XN_In(api,ip,e1i) + (1.0_rk - apw) * (1.0_rk - pD_In(api+1, ip, e1i)) * XN_In(api+1,ip,e1i)
			xiRVal = xiRVal + xiProb(j) * iTrans(i, ip) * (1.0_rk - xiGrid(j)) * (DVal + NVal)
		end do
	end do 
	xiRVal = xiRval / (1.0_rk + rf)
	xiR(ai,i,e1i) = xiRVal 
	xiTildeR(ai,i,e1i) = xiRVal * (-aGrid(ai))
	xijVec = (/ xi7_In(ai, e1i, e2i, e3i), xi13_In(ai, e1i, e2i, e3i), xiRVal /)	! vector of recovery rates by default type 
	TXD(ai, i, e1i) = dot_product(xijVec, pDjVec)
end do
!$omp end parallel do 
end do
end do 

! account for all debt actions
do e1i = 1, ne1
do i = 1,ni
!$omp parallel do private(LoanSize, TqVal, pDVal)
do ai = 1,na
do api = 1,naNeg
	LoanSize = LoanSizeInt(ai, api)
	if (LoanSize .gt. 0.0_rk) then 
		TqVal = -LoanFixedCost / LoanSize		! if its a net loan, account for fixed cost before default
	else
		TqVal = 0.0_rk 
	end if 
	
	! compute weighted expectation across future states 
	do ip = 1,ni 
		pDVal = pD_In(api,ip,e1i)
		TqVal = TqVal + iTrans(i, ip) * (pDVal * TXD(api, ip, e1i) + (1.0_rk - pDVal) * TXN(api, ip, e1i)) / (1.0_rk + rf + IntermediationCost)
	end do 
	Tq(api, ai, i, e1i) = TqVal
end do
end do
!$omp end parallel do
end do
end do 

end subroutine UpdatePriceFunction











!======================================================!
!                                                      !
! COMPUTE CONSUMPTION AND SUNNY DAY ADJUSTMENT WEIGHTS !
!                                                      !
!======================================================!
! NB: this is divided into 2 subroutines: 
! - one for all savings actions, which does not depend on flag status or the current price vector
! - one for all borrowing actions, which is only for no flag and does not depend on the current price vector

! SUBROUTINE 1: SAVING 
subroutine ComputeSunnyDay_Save(cVec, dVec, apiMaxPos, CantSave, yVal, aVal)
real(rk), intent(in) :: yVal, aVal
integer(ik), intent(out) :: apiMaxPos, CantSave
real(rk), intent(out), dimension(naPos) :: cVec, dVec
real(rk) :: apMaxVal

apMaxVal = (yVal + (kappa + (1.0_rk - phi) * qBar) * aVal - cMin) / qBar	! maximum feasible savings
cVec = 0.0_rk 
dVec = 0.0_rk 

if (apMaxVal .gt. 0.0_rk) then 
	CantSave = 0_ik 
	if (apMaxVal .gt. aGridPos(naPos)) then 
		apiMaxPos = naPos
	else 
		apiMaxPos = gridlookup(naPos, aGridPos, apMaxVal)
	end if
	cVec(1:apiMaxPos) = yVal + (kappa + (1.0_rk - phi) * qBar) * aVal - qBar * aGridPos(1:apiMaxPos)

	! will have to adjust for the no flag case at the 0 point
	dVec = 0.0_rk
	dVec(1) = (cVec(1) - cVec(2)) / 2.0_rk 										!	c(a0) - (c(a0) + c(a1)) / 2; note cVec(aZeroInd) = CashOnHand by definition
	dVec(2:apiMaxPos-1) = (cVec(1:apiMaxPos-2) - cVec(3:apiMaxPos)) / 2.0_rk	! 	(c(ai-1) + (c(ai)) / 2 - (c(ai) + c(ai+1)) / 2
	dVec(apiMaxPos) = (cVec(apiMaxPos) + cVec(apiMaxPos-1)) / 2.0_rk - cMin		! 	remaining span of budget set 
else
	CantSave = 1_ik 	! this household cannot save, and the feasible set will have to be adjusted.
	apiMaxPos = 1_ik
end if 
end subroutine ComputeSunnyDay_Save






! SUBROUTINE 2: BORROW 
subroutine ComputeSunnyDay_Borrow(cVec, dVec, apiMin, MustDefault, CantSave, apiMax, c0, c1, yVal, aVal, qVec, ai)
real(rk), intent(in) :: yVal, aVal, c0, c1	! c0 and c1 are the consumption at a' = 0 and the first grid point to the right of that, which are needed to compute the weights on the actions around 0 here.
real(rk), intent(in), dimension(naNeg) :: qVec
integer(ik), intent(in) :: CantSave, ai

integer(ik), intent(inout) :: apiMax

integer(ik), intent(out) :: apiMin, MustDefault
real(rk), intent(out), dimension(naNeg) :: cVec
real(rk), intent(out), dimension(naNeg+1) :: dVec	! has to be a bit longer in order to include 0

real(rk) :: apMaxVal
real(rk), dimension(naNeg) :: qhatVec
real(rk), dimension(naNeg+2) :: cVecTmp
integer(ik) :: api, apiLast, apiNext

! consumption associated with all debt actions 
if (aVal .lt. 0.0_rk) then
	qhatVec = qVec		! if in debt, prices governed totally by choice 
else
	qhatVec = qBar		! if not in debt , can always sell back at risk-free price
end if
cVec = yVal + (kappa + (1.0_rk - phi) * qhatVec) * aVal - qVec * aGridNeg

where (cVec .lt. cMin)
	cVec = 0.0_rk 
end where 


dVec = 0.0_rk
if (maxval(cVec) .lt. cMin) then 		! in this case, there are no feasible actions and the borrower must default (by definition, no feasible savings actions in this case)
	MustDefault = 1_ik
	cVec = 0.0_rk
else 
	MustDefault = 0_ik
	cVecTmp(1:naNeg) = cVec 
	cVecTmp(naNeg+1) = c0 
	cVecTmp(naNeg+2) = c1
	apiMin = maxloc(cVec, 1)	! peak of the Laffer Curve, make this the minimum feasible action
	
	if (CantSave .eq. 1_ik) then 	! adjust upper bound of feasible a' set to reflect the fact that borrowing must occur.
		apiMax = apiMin 
		do while (cVecTmp(apiMax+1) .gt. cMin)
			apiMax = apiMax + 1_ik
		end do 
	end if 	! otherwise apiMax remains the same from the "save" routine

	! work backwards through the debt actions to compute the weights correctly. 
	api = aZeroInd
	apiLast = aZeroInd + 1_ik
	apiNext = aZeroInd - 1_ik
	do while (apiNext .ge. apiMin)
		if (cVecTmp(apiNext) .gt. cVecTmp(api)) then 	! "normal" case: going deeper into debt yields more current consumption
			dVec(api) =  (cVecTmp(apiNext) - cVecTmp(apiLast)) / 2.0_rk	
			apiLast = api
			api = apiNext
		else 									! "pathological" case: going deeper into debt yields less current consumption, can't be optimal.
			dVec(apiNext) = 0.0_rk
		end if 
		apiNext = apiNext - 1_ik
	end do
	dVec(apiMin) = (cVecTmp(apiMin) - cVecTmp(api)) / 2.0_rk
	dVec(api) = (cVecTmp(apiMin) - cVecTmp(apiLast)) / 2.0_rk
end if 

end subroutine ComputeSunnyDay_Borrow






! COMPUTE CONSUMPTION, DECISION WEIGHTS, AND MAXIMAL INDEX FOR ALL SAVINGS ACTIONS AT ALL POINTS IN THE STATE SPACE
! this is called outside of the partial equilibrium loop since it does not depend on the loan price vector,
! 	but then must be "stitched together" to form the full action- / state-specific matrices one the loan price vector is fixed.
subroutine ConsumptionAndWeights_Save()

integer(ik) :: ai, i, e1i, li
real(rk) :: yVal 
integer(ik), dimension(na,ni,ne1,nl) :: apiMaxPosMat
real(rk), dimension(naPos,na,ni,ne1,nl) :: cSaveMat, dSaveMat

do li = 1,nl
do e1i = 1,ne1
do i = 1,ni
	yVal = yMat(e1i, e2Map(i), e3Map(i)) * (1.0_rk - lambdaGrid(li))
!$omp parallel do
do ai = 1,na
	call ComputeSunnyDay_Save(cSaveMat(:,ai,i,e1i,li), dSaveMat(:,ai,i,e1i,li), apiMaxPosMat(ai,i,e1i,li), CantSave0(ai,i,e1i), yVal, aGrid(ai))
end do
!$omp end parallel do
end do
end do
end do

! ASSIGN TO THE TYPE-SPECIFIC OBJECTS:
! consumption
CA0(aZeroInd:na,:,:,:) = cSaveMat(:,:,:,:,1)	! all current a, lambda = 0 for no flag, leave room for debt actions 
CA7 = cSaveMat(:,aZeroInd:na,:,:,1)						! only positive current a, lambda = 0 for Ch 7
CA13 = cSaveMat(:,aZeroInd:na,:,:,:)					! only positive current a, all lambda for Ch 13

! action weights 
dA0(aZeroInd+1:na,:,:,:) = dSaveMat(2:naPos,:,:,:,1)	! have to allow a' = 0 weight to be written over based on available consumption in debt
dA7 = dSaveMat(:,aZeroInd:na,:,:,1)	
dA13 = dSaveMat(:,aZeroInd:na,:,:,:)

! maximal feasible a' index: note that the "scale" of the maximal index reflects the feasible choice range given the flag status (naPos for Ch 7 or 13, na for no flag)
apiMax0 = aiFull(apiMaxPosMat(:,:,:,1))
apiMax7 = apiMaxPosMat(aZeroInd:na,:,:,1)
apiMax13 = apiMaxPosMat(aZeroInd:na,:,:,:)

end subroutine ConsumptionAndWeights_Save







! this is called inside the partial equilibrium loop because it depends on loan prices.
! NB: it also depends on the first few values of consumption computed in the saving version
subroutine ConsumptionAndWeights_Borrow(q_In)
real(rk), intent(in), dimension(naNeg, na, ni, ne1) :: q_In		! current price schedule
integer(ik) :: ai, i, e1i

do e1i = 1,ne1
do i = 1,ni 
do ai = 1,na 
	call ComputeSunnyDay_Borrow(CA0(1:naNeg,ai,i,e1i), dA0(1:aZeroInd,ai,i,e1i), apiMin0(ai,i,e1i), MustDefault0(ai,i,e1i), &	! outputs
		& CantSave0(ai,i,e1i), apiMax0(ai,i,e1i), &
		& CA0(aZeroInd,ai,i,e1i), CA0(aZeroInd+1,ai,i,e1i), yMat(e1i, e2Map(i), e3Map(i)), aGrid(ai), q_In(:,ai,i,e1i), ai)			! inputs 
end do
end do
end do 

end subroutine ConsumptionAndWeights_Borrow









!=========================!
!                         !
! STATIONARY DISTRIBUTION !
!                         !
!=========================!
subroutine ComputeAggregateCapital(KVal, m0, m7, m13)
real(rk), intent(out) :: KVal
real(rk), intent(in), dimension(na, ni, ne1) :: m0
real(rk), intent(in), dimension(naPos, ni, ne1) :: m7
real(rk), intent(in), dimension(naPos, ni, ne1, nl) :: m13
integer(ik) :: ai, aiPos 
real(rk) :: aVal 

KVal = 0.0_rk
do ai = 1,na
	aVal = aGrid(ai)
	KVal = KVal + sum(mu0(ai,:,:)) * aVal
	aiPos = aiPositive(ai)
	if (aiPos .ge. 1_ik) then 
		KVal = KVal + sum(mu7(aiPos,:,:)) * aVal
		KVal = KVal + sum(mu13(aiPos,:,:,:)) * aVal
	end if 
end do 
end subroutine ComputeAggregateCapital




! SOLVE FOR STATIONARY DISTRIBUTION
subroutine sdSolve(m0, m7, m13, pA0_In, pA7_In, pA13_In, pD_In, pD7_In, pD13_In, pDR_In, pA0_GoDQ_In)
real(rk), intent(in), dimension(na, na, ni, ne1) :: pA0_In
real(rk), intent(in), dimension(naPos, naPos, ni, ne1) :: pA7_In
real(rk), intent(in), dimension(naPos, naPos, ni, ne1, nl) :: pA13_In

real(rk), intent(in), dimension(naNeg, ni, ne1) :: pD_In, pD7_In, pD13_In, pDR_In
real(rk), intent(in), dimension(naNeg+1, naNeg) :: pA0_GoDQ_In

real(rk), intent(inout), dimension(na, ni, ne1) :: m0
real(rk), intent(inout), dimension(naPos, ni, ne1) :: m7
real(rk), intent(inout), dimension(naPos, ni, ne1, nl) :: m13

integer(ik) :: mIter
real(rk)    :: mDist, m0Dist, m7Dist, m13Dist, mSum, TmSum
real(rk), dimension(na, ni, ne1) :: Tm0, Tm0_Ch7, Tm0_Ch13
real(rk), dimension(naPos, ni, ne1) :: Tm7, Tm7_Ch7
real(rk), dimension(naPos, ni, ne1, nl) :: Tm13, Tm13_Ch13

! MAIN CALCULATIONS 
mDist = 1.0_rk
mIter = 0_ik
144 format(i10,f10.6, 10f10.4)
do while (mDist .gt. muTol .and. mIter .lt. muMaxIters)
    mIter = mIter + 1_ik
	
	! compute how the distribution updates for each current flag status 
	call sdSingleUpdate_0(Tm0, Tm7, Tm13, m0, pA0_In, pD_In, pD7_In, pD13_In, pDR_In, pA0_GoDQ_In, 1, na)	! no flag 
	call sdSingleUpdate_Ch7(Tm0_Ch7, Tm7_Ch7, m7, pA7_In, 1, na)											! Chapter 7
	call sdSingleUpdate_Ch13(Tm0_Ch13, Tm13_Ch13, m13, pA13_In, 1, na)										! Chapter 13
	
	! accumulate across subroutines 
	Tm0 = Tm0 + Tm0_Ch7 + Tm0_Ch13 
	Tm7 = Tm7 + Tm7_Ch7
	Tm13 = Tm13 + Tm13_Ch13
	
	! assess convergence and update
    m0Dist = maxval(abs(m0 - Tm0))
    m7Dist = maxval(abs(m7 - Tm7))
    m13Dist = maxval(abs(m13 - Tm13))
	mDist = maxval((/ m0Dist, m7Dist, m13Dist /))
	mSum = sum(m0) + sum(m7) + sum(m13)
	TmSum = sum(Tm0) + sum(Tm7) + sum(Tm13)
    if (mod(mIter, 20_ik) .eq. 0_ik .or. mIter .le. 5_ik) then
        print 144, mIter, mDist, mSum, TmSum, sum(Tm0), sum(Tm7), sum(Tm13)
    end if
    m0 = Tm0
    m7 = Tm7 
    m13 = Tm13
	
end do      ! end of stationary distribution iteration loop

print*, ' '
print 144, mIter, mDist, sum(m0) + sum(m7) + sum(m13)
print*, ' '

end subroutine sdSolve






! ITERATE ON T* OPERATOR TO SOLVE FOR STATIONARY DISTRIBUTION
subroutine sdSingleUpdate_Ch7(Tm0, Tm7, m7_In, pA7_In, ai1, ai2)
! this subroutine only accounts for the behavior and transitions of those CURRENTLY with the Ch 7 flag
integer(ik), intent(in) :: ai1, ai2
real(rk), intent(out), dimension(na, ni, ne1) :: Tm0
real(rk), intent(out), dimension(naPos, ni, ne1) :: Tm7
real(rk), intent(in), dimension(naPos, ni, ne1) :: m7_In
real(rk), intent(in), dimension(naPos, naPos, ni, ne1) :: pA7_In

integer(ik) :: aiPos, i, e1i, ai, ip
real(rk) :: mVal
real(rk), dimension(naPos) :: probVec

Tm0 = 0.0_rk 	! this will have to be added back to the ones implied by the others 
Tm7 = 0.0_rk 	! this will have to reflect new flows in from flag = 0 as well 
do e1i = 1, ne1
do i = 1, ni
do aiPos = 1,naPos
	mVal = m7_In(aiPos, i, e1i)
	if (mVal .gt. 0.0_rk) then
		do ip = 1,ni
			probVec = iTrans(i, ip) * mVal * pA7_In(:,aiPos,i,e1i)
			Tm7(:,ip,e1i) = Tm7(:,ip,e1i) + (1.0_rk - theta_7) * probVec
			Tm0(aZeroInd:na,ip,e1i) = Tm0(aZeroInd:na,ip,e1i) + theta_7 * probVec
		end do
	end if
end do  
end do  
end do  

end subroutine sdSingleUpdate_Ch7




subroutine sdSingleUpdate_Ch13(Tm0, Tm13, m13_In, pA13_In, ai1, ai2)
! this subroutine only accounts for the behavior and transitions of those CURRENTLY with the Ch 13 flag
integer(ik), intent(in) :: ai1, ai2
real(rk), intent(out), dimension(na, ni, ne1) :: Tm0
real(rk), intent(out), dimension(naPos, ni, ne1, nl) :: Tm13
real(rk), intent(in), dimension(naPos, ni, ne1, nl) :: m13_In
real(rk), intent(in), dimension(naPos, naPos, ni, ne1, nl) :: pA13_In
integer(ik) :: aiPos, i, e1i, ai, li, ip
real(rk) :: mVal
real(rk), dimension(naPos) :: probVec

Tm0 = 0.0_rk
Tm13 = 0.0_rk 
do li = 1,nl
do e1i = 1, ne1
do i = 1, ni
do aiPos = 1,naPos
	mVal = m13_In(aiPos, i, e1i, li)
	if (mVal .gt. 0.0_rk) then
		do ip = 1,ni
			probVec = iTrans(i, ip) * mVal * pA13_In(:,aiPos,i,e1i,li)
			Tm13(:,ip,e1i,li) = Tm13(:,ip,e1i,li) + (1.0_rk - theta_13) * probVec
			Tm0(aZeroInd:na,ip,e1i) = Tm0(aZeroInd:na,ip,e1i) + theta_13 * probVec
		end do
	end if
end do  ! ai
end do  ! i
end do  ! e1i
end do 	! li

end subroutine sdSingleUpdate_Ch13




subroutine sdSingleUpdate_0(Tm0, Tm7, Tm13, m0_In, pA0_In, pD_In, pD7_In, pD13_In, pDR_In, pA0_GoDQ_In, ai1, ai2)
integer(ik), intent(in) :: ai1, ai2
real(rk), intent(out), dimension(na, ni, ne1) :: Tm0
real(rk), intent(out), dimension(naPos, ni, ne1) :: Tm7
real(rk), intent(out), dimension(naPos, ni, ne1, nl) :: Tm13
real(rk), intent(in), dimension(na, ni, ne1) :: m0_In
real(rk), intent(in), dimension(na, na, ni, ne1) :: pA0_In
real(rk), intent(in), dimension(naNeg, ni, ne1) :: pD_In, pD7_In, pD13_In, pDR_In
real(rk), intent(in), dimension(naNeg+1, naNeg) :: pA0_GoDQ_In
integer(ik) :: ai, i, e1i, api, ip, lpi, e2i, e3i
real(rk) :: mVal, prob, apw, pDVal, lpw, pD7Val, pD13Val, pDRVal
real(rk), dimension(na) :: NashProbVec 

Tm0 = 0.0_rk
Tm7 = 0.0_rk
Tm13 = 0.0_rk
do e1i = 1, ne1
do i = 1, ni
	e2i = e2Map(i)
	e3i = e3Map(i)
do ai = ai1, ai2
	mVal = m0_In(ai, i, e1i)
	if (mVal .gt. 0.0_rk) then
		if (ai .lt. aZeroInd) then 			! default is possible, have to adjust for this 
			pDVal = pD_In(ai, i, e1i)		! probability of default of any type
			pD7Val = pD7_In(ai, i, e1i)		! Ch 7 prob 
			pD13Val = pD13_In(ai, i, e1i)	! Ch 13 prob 
			lpi = lpiMat(ai, e1i, e2i, e3i)	! index and weight of the garnishment fraction going forward in Ch 13
			lpw = lpwMat(ai, e1i, e2i, e3i)
			pDRVal = pDR_In(ai, i, e1i)		! restructuring prob 
			
			do ip = 1,ni   
				prob = iTrans(i, ip) * mVal  
				Tm0(:,ip,e1i) = Tm0(:,ip,e1i) + prob * (1.0_rk - pDVal) * pA0_In(:,ai,i,e1i)				! no default
				Tm7(1,ip,e1i) = Tm7(1,ip,e1i) + prob * pDVal * pD7Val 										! Chapter 7
				Tm13(1,ip,e1i,lpi) = Tm13(1,ip,e1i,lpi) + prob * pDVal * pD13Val * lpw
				Tm13(1,ip,e1i,lpi+1) = Tm13(1,ip,e1i,lpi+1) + prob * pDVal * pD13Val * (1.0_rk - lpw)		! Chapter 13
				Tm0(1:naNeg+1,ip,e1i) = Tm0(1:naNeg+1,ip,e1i) + prob * pDVal * pDRVal * pA0_GoDQ_In(:,ai)	! restructuring / delinquency
			end do
		else 						! default is not possible, just borrow / save decisions 
			do ip = 1,ni   
				Tm0(:,ip,e1i) = Tm0(:,ip,e1i) + iTrans(i, ip) * mVal * pA0_In(:,ai,i,e1i)
			end do 
		end if 
	end if  ! end of default / no default conditional for people with a Flag.
end do  ! ai
end do  ! i
end do  ! e1i
end subroutine sdSingleUpdate_0







!==============================================!
!                                              !
! COMPUTE STEADY STATE MOMENTS AND SAVE OUTPUT !
!                                              !
!==============================================!
! COMPUTE TARGET MOMENTS 
subroutine ComputeSSMoments()
integer(ik) :: ai, i, e1i, ip, api, e2i, e3i, aiPos
real(rk) :: muVal, aVal, apVal, yVal, pDVal, pD7Val, pD13Val, pDRVal, pBKVal, AverageSpreadSquared, TotalDebt, pAVal, TotalLoans, SpreadDenom
real(rk) :: YTM_riskfree, YTMVal, SpreadVal, CORVal, RRVal, LoanVal, aMean, TotalDollarsWithDefault, TotalDollarsLost, TotalDollarsWithBK, TotalDollarsLostBK,TotalDollarsDebt, xi7Val, xi13Val, xiRVal 
real(rk), dimension(2) :: RegMoments 

real(rk), dimension(naNeg, ni, ne1) :: muDef, muCh7, muCh13, muR
real(rk) :: aBarCh13, aBarCh7, aBarDef, aBarR, yBarCh13, yBarCh7, yBarDef, yBarR, rBarCh13, rBarCh7, rBarDef, rBarR
real(rk), dimension(4,4) :: X 

! INITIALIZATIONS 
DebtToIncome = 0.0_rk    			! avg debt to labor income (savers have 0)
DebtToIncome_ap = 0.0_rk
DefaultRate = 0.0_rk    			! share of population defaulting within a period
ChargeoffRate = 0.0_rk    			! share of total dollars defaulted on within a period
RecoveryRate = 0.0_rk    			! share of total dollars defaulted on within a period
Ch7Rate = 0.0_rk
Ch13Rate = 0.0_rk
RestructuringRate = 0.0_rk 

SpreadDenom = 0.0_rk
TotalDebt = 0.0_rk    				! economywide total debt
TotalLoans = 0.0_rk 
AverageSpread = 0.0_rk     			! avg interest rate on all borrowing
AverageSpreadSquared = 0.0_rk 		! avg interest rate  squared on all borrowing
TotalDollarsWithDefault = 0.0_rk
TotalDollarsLost = 0.0_rk
TotalDollarsWithBK = 0.0_rk
TotalDollarsLostBK = 0.0_rk
TotalDollarsDebt = 0.0_rk

YTM_riskfree = kappa - phi			! = r, YTM on risk-free bond (given normalization of qrf to 1)

FractionInDebt = sum(mu0(1:naNeg,:,:))

! LOOP OVER ALL STATES AND CHOICES
do e1i = 1, ne1
do i = 1, ni
	e2i = e2Map(i)
	e3i = e3Map(i)
	yVal = yMat(e1i, e2i, e3i)
do ai = 1,na
	aVal = aGrid(ai)
	
    muVal = mu0(ai,i,e1i)
	if (muVal .gt. 0.0_rk) then 
		if (ai .lt. aZeroInd) then 
			pDVal = pD(ai,i,e1i)
			
			pD7Val = pD7(ai,i,e1i)
			xi7Val = xi7(ai, e1i, e2i, e3i)
			
			pD13Val = pD13(ai,i,e1i)
			xi13Val = xi13(ai,e1i,e2i,e3i)
			
			pDRVal = pDR(ai,i,e1i)
			xiRVal = xiR(ai, i, e1i)
			
			pBKVal = pD7Val + pD13Val 
			
			DefaultRate = DefaultRate + muVal * pDVal	
			Ch7Rate = Ch7Rate + muVal * pDVal * pD7Val	
			Ch13Rate = Ch13Rate + muVal * pDVal * pD13Val
			RestructuringRate = RestructuringRate + muVal * pDVal * pDRVal
			if (ExpDefProb(ai,i,e1i,nAhead) .le. pDMax) then 
				DebtToIncome = DebtToIncome + muVal * (-aVal / yVal)
			end if 
			
			TotalDollarsLost = TotalDollarsLost + muVal * (-aVal) * pDVal * (pD7Val * (1.0_rk - xi7Val) + pD13Val * (1.0_rk - xi13Val) + pDRVal * (1.0_rk - xiRval))
			TotalDollarsWithDefault = TotalDollarsWithDefault + muVal * (-aVal) * pDVal
			
			TotalDollarsLostBK = TotalDollarsLostBK + muVal * (-aVal) * pDVal * (pD7Val * (1.0_rk - xi7Val) + pD13Val * (1.0_rk - xi13Val))
			TotalDollarsWithBK = TotalDollarsWithBK + muVal * (-aVal) * pDVal * (pD7Val + pD13Val)
			
			TotalDollarsDebt = TotalDollarsDebt + muVal * (-aVal)
		end if 
		
		! loop over all debt choices to accumulate moments as necessary
		do api = 1,naNeg
			apVal = aGrid(api)
			if (ai .lt. aZeroInd) then 
				pAVal = (1.0_rk - pDVal) * pA0(api,ai,i,e1i)
			else
				pAVal = pA0(api,ai,i,e1i)
			end if 
			if (pAVal .gt. 0.0_rk) then 
				LoanVal = LoanSizeInt(ai, api)
				spreadVal = (1.0_rk + spreads(api,ai,i,e1i)) ** 4.0_rk - 1.0_rk
				
				TotalLoans = TotalLoans + LoanVal * muVal * pAVal 
				TotalDebt = TotalDebt - muVal * apVal * pAVal 
				
				if (ExpDefProb(ai,i,e1i,nAhead) .le. pDMax) then 
					AverageSpread = AverageSpread - spreadVal * muVal * pAval * apVal
					AverageSpreadSquared = AverageSpreadSquared - (spreadVal ** 2.0_rk) * muVal * pAval * apVal
					SpreadDenom = SpreadDenom -  muVal * pAval * apVal
					DebtToIncome_ap = DebtToIncome_ap - muVal * pAVal * (apVal / yVal)
				end if 
			end if 
		end do 
	end if 
end do
end do
end do

print*, 'DTI ap',DebtToIncome_ap

! scale as needed for each moment
RecoveryRate = (TotalDollarsWithBK - TotalDollarsLostBK) / TotalDollarsWithBK
ChargeoffRate = TotalDollarsLost / TotalDollarsDebt
BankruptcyRate = Ch7Rate + Ch13Rate
Ch7Share = Ch7Rate / BankruptcyRate

if (SpreadDenom .gt. 0.0_rk) then 
	AverageSpread = AverageSpread / SpreadDenom
	AverageSpreadSquared = AverageSpreadSquared / SpreadDenom
	SDSpread = sqrt(AverageSpreadSquared - AverageSpread ** 2.0_rk) 
else 
	AverageSpread = 0.0_rk
	SDSpread = 0.0_rk
end if 

aMean = 0.0_rk 
do ai = 1,na
	aMean = aMean + aGrid(ai) * sum(mu0(ai,:,:))
	if (ai .ge. aZeroInd) then 
		aiPos = aiPositive(ai)
		aMean = aMean + aGrid(ai) * (sum(mu7(aiPos,:,:)) + sum(mu13(aiPos,:,:,:)))
	end if 
end do 
print*, ' '
print*, 'mean A:', aMean
print*, ' '



!--------------------------------------------------------!
! default and recovery statistics (response to Cristina) !
!--------------------------------------------------------!
! partial distributions 
muDef = mu0(1:naNeg,:,:) * pD 			! distribution of defaulters 
muCh7 = mu0(1:naNeg,:,:) * pD * pD7			! distribution of Ch 7 bankrutpts
muCh13 = mu0(1:naNeg,:,:) * pD * pD13		! distribution of Ch 13 bankrutpts
muR = mu0(1:naNeg,:,:) * pD * pDR			! distribution of Ch delinquents

! average assets
aBarDef = 0.0_rk 
aBarCh7 = 0.0_rk 
aBarCh13 = 0.0_rk 
aBarR = 0.0_rk  
do ai = 1,naNeg 
	aVal = -aGrid(ai)
	aBarDef = aBarDef + aVal * sum(muDef(ai,:,:)) / sum(muDef)
	aBarCh7 = aBarCh7 + aVal * sum(muCh7(ai,:,:)) / sum(muCh7)
	aBarCh13 = aBarCh13 + aVal * sum(muCh13(ai,:,:)) / sum(muCh13)
	aBarR = aBarR + aVal * sum(muR(ai,:,:)) / sum(muR)
end do 

! average income 
yBarDef = 0.0_rk 
yBarCh7 = 0.0_rk 
yBarCh13 = 0.0_rk 
yBarR = 0.0_rk  
do e1i = 1,ne1
do i = 1,ni 
	e2i = e2Map(i)
	e3i = e3Map(i)
	yVal = yMat(e1i, e2i, e3i)
	
	yBarDef = yBarDef + yVal * sum(muDef(:,i,e1i)) / sum(muDef)
	yBarCh7 = yBarCh7 + yVal * sum(muCh7(:,i,e1i)) / sum(muCh7)
	yBarCh13 = yBarCh13 + yVal * sum(muCh13(:,i,e1i)) / sum(muCh13)
	yBarR = yBarR + yVal * sum(muR(:,i,e1i)) / sum(muR)
end do 
end do 


! average recovery 
rBarDef = 0.0_rk 
rBarCh7 = 0.0_rk 
rBarCh13 = 0.0_rk 
rBarR = 0.0_rk  
do e1i = 1,ne1
do i = 1,ni 
	e2i = e2Map(i)
	e3i = e3Map(i)
do ai = 1,naNeg 
	xi7Val = xi7(ai, e1i, e2i, e3i)
	xi13Val = xi13(ai, e1i, e2i, e3i)
	xiRVal = xiR(ai, i, e1i)
	
	rBarCh7 = rBarCh7 + xi7(ai, e1i, e2i, e3i) * muCh7(ai,i,e1i) / sum(muCh7)
	rBarCh13 = rBarCh13 + xi13(ai,e1i,e2i,e3i) * muCh13(ai,i,e1i) / sum(muCh13)
	rBarR = rBarR + xiR(ai, i, e1i) * muR(ai,i,e1i) / sum(muR)
	
	rBarDef = rBarDef + (pD7(ai,i,e1i) * xi7Val + pD13(ai,i,e1i) * xi13Val + pDR(ai, i, e1i) * xiRVal) * muDef(ai,i,e1i) / sum(muDef)
end do 
end do 
end do 

! compile
X(1,:) = (/ 100.0_rk * ((1.0_rk + DefaultRate) ** 4.0_rk - 1.0_rk), aBarDef, yBarDef, 100.0_rk * rBarDef /)
X(2,:) = (/ 100.0_rk * ((1.0_rk + Ch7Rate) ** 4.0_rk - 1.0_rk), aBarCh7, yBarCh7, 100.0_rk * rBarCh7 /)
X(3,:) = (/ 100.0_rk * ((1.0_rk + Ch13Rate) ** 4.0_rk - 1.0_rk), aBarCh13, yBarCh13, 100.0_rk * rBarCh13 /)
X(4,:) = (/ 100.0_rk * ((1.0_rk + RestructuringRate) ** 4.0_rk - 1.0_rk), aBarR, yBarR, 100.0_rk * rBarR /)

! display 
315 format(a15,5f10.2)
print*, ' '
print 315, 'any default', X(1,:)
print 315, 'Ch 7 BK', X(2,:)
print 315, 'Ch 13 BK', X(3,:)
print 315, 'restructuring', X(4,:)
print*, ' '

open (unit = 1, file = trim(OutDir)//trim('Table8'), status = 'replace')
do i = 1,4
	write (unit = 1, fmt = "(10f8.2)") X(i,:)
end do
close (unit = 1)




!---------------------------------------------------------!
! compute the spreads by default probability relationship !
!---------------------------------------------------------!
print*,' allocating partition...'
call linspace(pDMin, pDMax, nPartition, pDPartition)	! evenly spaced partition of default probabilities
pDPartitionPlot(1) = pDPartition(1) / 2.0_rk 			! adjust to find midpoints 
do i = 2,nPartition
	pDPartitionPlot(i) = (pDPartition(i-1) + pDPartition(i)) / 2.0_rk 			! adjust to find midpoints 
end do

print*, 'computing bin-specific spreads...'
if (nAhead .eq. 1_ik) then 
	call AverageSpreadsPartition(RegMoments, ExpDefProb1)
else 
	call AverageSpreadsPartition(RegMoments, ExpDefProb(:,:,:,nAhead))
end if

SpreadIntercept = RegMoments(1)
SpreadSlope = RegMoments(2)
 
end subroutine ComputeSSMoments









!========================================================================!
!                                                                        !
! MAPPING SPREADS BY DEFAULT PROBABILITY INTO THE EMPIRICAL RELATIONSHIP !
!                                                                        !
!========================================================================!

! CONDITIONAL: COMPUTE THE ONE-PERIOD AHEAD DEFAULT PROBABILITY GIVEN THE A' CHOICE, THE CURRENT STATE, AND KNOWLEDGE OF NO CURRENT DEFAULT 
subroutine Conditional_OnePeriodAheadDefaultProbability(delta_1)
real(rk), intent(out), dimension(naNeg, ni, ne1)  :: delta_1 
integer(ik) :: e1i, i, api
do e1i = 1,ne1 
do i = 1,ni 
	!$omp parallel do
	do api = 1,naNeg
		delta_1(api, i, e1i) = dot_product(iTrans(i,:), pD(api,:,e1i))
	end do
	!$omp end parallel do
end do 
end do 
end subroutine Conditional_OnePeriodAheadDefaultProbability

! ITERATIVELY COMPUTE THE N-PERIOD AHEAD CONDITIONAL DEFAULT PROBABILITY GIVEN THE N-1 DEFAULT PROBABILITY AND THE COMBINED EXOGENOUS / ENDOGENOUS TRANSITIONS 
subroutine Conditional_NPeriodAheadDefaultProbability(delta_i, delta_im1, delta_1)
real(rk), intent(in), dimension(naNeg, ni, ne1)  :: delta_im1, delta_1
real(rk), intent(out), dimension(naNeg, ni, ne1)  :: delta_i
real(rk) :: tmp 
integer(ik) :: e1i, i, api, appi 

do e1i = 1,ne1
do i = 1,ni 
!$omp parallel do private(tmp)
do api = 1,na 
	tmp = delta_1(api, i, e1i)
	do appi = 1,na
		tmp = tmp + sum(iTrans(i,:) * delta_im1(appi, :, e1i) * pA0(appi, api, :, e1i) * (1.0_rk - pD(api,:,e1i)))
	end do
	delta_i(api, i, e1i) = tmp 
end do
!$omp end parallel do 
end do
end do 
end subroutine Conditional_NPeriodAheadDefaultProbability



! UNCONDITIONAL: COMPUTE THE ONE-PERIOD AHEAD DEFAULT PROBABILITY GIVEN ONLY CURRENT STATE AND KNOWLEDGE OF NO CURRENT DEFAULT
subroutine OnePeriodAheadDefaultProbability(delta_1)
real(rk), intent(out), dimension(na, ni, ne1)  :: delta_1 
integer(ik) :: e1i, i, ai, api 
real(rk) :: tmp 

do e1i = 1,ne1 
do i = 1,ni 
!$omp parallel do private(tmp)
do ai = 1,na
	tmp = 0.0_rk 
	do api = 1,naNeg
		tmp = tmp + pA0(api, ai, i, e1i) * dot_product(iTrans(i,:), pD(api,:,e1i))
	end do
	delta_1(ai, i, e1i) = tmp
end do 
!$omp end parallel do
end do 
end do 

end subroutine OnePeriodAheadDefaultProbability


! ITERATIVELY COMPUTE THE N-PERIOD AHEAD DEFAULT PROBABILITY GIVEN THE N-1 DEFAULT PROBABILITY AND THE COMBINED EXOGENOUS / ENDOGENOUS TRANSITIONS 
subroutine NPeriodAheadDefaultProbability(delta_i, delta_im1, delta_1)
real(rk), intent(in), dimension(na, ni, ne1)  :: delta_im1, delta_1
real(rk), intent(out), dimension(na, ni, ne1)  :: delta_i
real(rk) :: tmp 
integer(ik) :: e1i, i, ai, api 

do e1i = 1,ne1
do i = 1,ni 
!$omp parallel do private(tmp)
do ai = 1,na 
	tmp = delta_1(ai, i, e1i)
	do api = 1,na
		tmp = tmp + pA0(api, ai, i, e1i) * sum(iTrans(i,:) * (1.0_rk - pD(api,:,e1i)) * delta_im1(api, :, e1i))
	end do
	delta_i(ai, i, e1i) = tmp
end do 
!$omp end parallel do 
end do
end do 


end subroutine NPeriodAheadDefaultProbability





! COMBINE SPREADS AND WEIGHTS TO COMPUTE SPREAD MOMENTS FOR EACH PARTITION INTERVAL 
subroutine AverageSpreadsPartition(SpreadProfile_Coefs, delta_n)
real(rk), intent(in), dimension(na, ni, ne1) :: delta_n	! the specified expected default probability metric
real(rk), intent(out), dimension(2) :: SpreadProfile_Coefs

integer(ik) :: e1i, i, api, j, ai
integer(ik), dimension(na, ni, ne1) :: jPartition
real(rk) :: mass, LoanSize, weight, minEpD, maxEpD, spreadVal, apVal, pAVal, mVal, pDVal
real(rk), dimension(nPartition) :: AvgBinSpreadSquared
real(rk), dimension(na, na, ni, ne1) :: MassChoiceState
real(rk), dimension(5) :: SpreadProfile_Stats

integer(ik) :: nReg 
integer(ik), allocatable, dimension(:) :: ones 
real(rk), allocatable, dimension(:) :: AvgBinSpreadReg, pDReg

! PARTITION THE CHOICE / STATE SPACE BASED ON DEFAULT PROBABILITY
MassPartition = 0.0_rk 
do e1i = 1,ne1
do i = 1,ni
do ai = 1,na 
	j = count(delta_n(ai, i, e1i) .gt. pDPartition) + 1_ik				! partition interval corresponding to default probability 
	jPartition(ai, i, e1i) = j
	mVal = mu0(ai, i, e1i)
	if (ai .lt. aZeroInd) then 
		pDVal = pD(ai, i, e1i)
	else
		pDVal = 0.0_rk 
	end if 
	MassPartition(j) = MassPartition(j) + (1.0_rk - pDVal) * dot_product(-aGrid(1:naNeg), pA0(1:naNeg, ai, i, e1i)) * mVal 
end do
end do
end do 

! at this point, the bin-specific weight on the action-state combo is MassChoiceState(a',a,i,e1) / MassPartition(jPartition(a',i,e1i))
AvgBinSpread = 0.0_rk 
AvgBinSpreadSquared = 0.0_rk 
do e1i = 1,ne1
do i = 1,ni
do ai = 1,na
	j = jPartition(ai, i, e1i)
	mVal = mu0(ai, i, e1i) 
	if (ai .lt. aZeroInd) then 
		pDVal = pD(ai, i, e1i)
	else
		pDVal = 0.0_rk 
	end if 
	do api = 1,naNeg
		weight = -aGrid(api) * (1.0_rk - pDVal) * pA0(api, ai, i, e1i) * mVal / MassPartition(j)
		if (weight .gt. 0.0_rk) then 
			spreadVal = 100.0_rk * ((1.0_rk + spreads(api,ai,i,e1i)) ** 4.0_rk - 1.0_rk)
			AvgBinSpread(j) = AvgBinSpread(j) + spreadVal * weight
			AvgBinSpreadSquared(j) = AvgBinSpreadSquared(j) + (spreadVal ** 2.0_rk) * weight
		end if 
	end do
end do
end do 
end do 

! adjust so that bins with no mass don't mess things up 
nReg = count(massPartition .gt. 0.0_rk)
allocate(ones(nReg), AvgBinSpreadReg(nReg), pDReg(nReg))
ones = 1_ik
i = 1_ik 
do j = 1,nPartition
	if (massPartition(j) .gt. 0.0_rk) then 
		AvgBinSpreadReg(i) = AvgBinSpread(j)
		pDReg(i) = 100.0_rk * pDPartition(j)
		i = i+1_ik
	end if 
end do 

! regress the profile 
SpreadProfile_Coefs = 0.0_rk
call invpolicy(nReg, 1_ik, 1_ik, ones, AvgBinSpreadReg, pDReg, SpreadProfile_Coefs, SpreadProfile_Stats)  

! report results on screen 
298 format(11a8)
299 format(i8, 10f8.3)
print*, ' '
print*, 'Average spreads by default probability '
print*, ' '
! print 298, 'j','min EpD','max EpD', 'mean s', 'fit', 'sd s', 'mass'
print 298, 'j','min EpD','max EpD', 'mean s', 'sd s', 'mass'
do j = 1,nPartition 
	StDevBinSpread(j) = sqrt(AvgBinSpreadSquared(j) - AvgBinSpread(j) ** 2.0_rk)
	if (j .eq. 1_ik) then
		minEpD = 0.0_rk 
		maxEpD = pDPartition(1)
	else
		minEpD = pDPartition(j-1)
		maxEpD = pDPartition(j)
	end if 
	print 299, j, minEpD, maxEpD, AvgBinSpread(j), StDevBinSpread(j), MassPartition(j)
end do 
print*, ' '
print*, 'Linear fit to profile'
300 format(a10,f10.3)
print 300, 'intercept', SpreadProfile_Coefs(1)
print 300, 'slope', SpreadProfile_Coefs(2)
print 300, 'R-squard', SpreadProfile_Stats(2)
print*, ' '

deallocate(AvgBinSpreadReg, pDReg, ones)


end subroutine AverageSpreadsPartition







!========!
!        !
! SAVING !
!        !
!========!
subroutine SaveSSOutput()

call savings(kappa,'kappa',OutDir)
call savings(phi,'phi',OutDir)

call savings(dble(na),'na',OutDir)
call savings(dble(naNeg),'naNeg',OutDir)
call savings(dble(naPos),'naPos',OutDir)
call savings(dble(ne1),'ne1',OutDir)
call savings(dble(ne2),'ne2',OutDir)
call savings(dble(ne3),'ne3',OutDir)
call savings(dble(ni),'ni',OutDir)
call savings(dble(nl),'nl',OutDir)
call savings(dble(nPartition),'nPartition',OutDir)

call savings(crra,'crra',OutDir)
call savings(ev_nu,'ev_nu',OutDir)
call savings(ev_zeta,'ev_zeta',OutDir)
call savings(ev_eta,'ev_eta',OutDir)
call savings(cMin,'cMin',OutDir)
call savings(alpha,'alpha',OutDir)
call savings(rf_ann,'rf_ann',OutDir)
call savings(macaulay,'macaulay',OutDir)
call savings(LoanFixedCost,'LoanFixedCost',OutDir)
call savings(stigma_7,'stigma_7',OutDir)
call savings(psi,'psi',OutDir)
call savings(xi_7,'xi_7',OutDir)
call savings(theta_7,'theta_7',OutDir)
call savings(stigma_13,'stigma_13',OutDir)
call savings(lambdaBar,'lambdaBar',OutDir)
call savings(xi_13,'xi_13',OutDir)
call savings(theta_13,'theta_13',OutDir)
call savings(stigma_R,'stigma_R',OutDir)
call savings(alpha_R,'alpha_R',OutDir)
call savings(qBar,'qBar',OutDir)
call savings(rf,'rf',OutDir)
call savings(wage,'wage',OutDir)

! value functions 
call saving3d(W0,'W0',OutDir)
call saving3d(WD0,'WD0',OutDir)
call saving3d(WN0,'WN0',OutDir)

! call saving3d(J7,'J7',OutDir)
! call saving3d(J13,'J13',OutDir)
! call saving3d(JR,'JR',OutDir)

call saving3d(W7,'W7',OutDir)
call saving4d(W13,'W13',OutDir)

! decision rules 
call saving4d(pA0,'pA0',OutDir)
call saving4d(pA7,'pA7',OutDir)
call saving5d(pA13,'pA13',OutDir)

! ! consumption
! call saving4d(CA0,'CA0',OutDir)
! call saving4d(CA7,'CA7',OutDir)
! call saving5d(CA13,'CA13',OutDir)

call saving3d(pd,'pD',OutDir)
call saving3d(pd7,'pD7',OutDir)
call saving3d(pd13,'pD13',OutDir)
call saving3d(pdR,'pDR',OutDir)

! distributions 
call saving3d(mu0,'mu0',OutDir)
call saving3d(mu7,'mu7',OutDir)
call saving4d(mu13,'mu13',OutDir)

! prices 
call saving4d(q,'q',OutDir)
call saving4d(spreads,'spreads',OutDir)
call saving4d(xi7,'xi7',OutDir)
call saving4d(xi13,'xi13',OutDir)
call saving3d(xiR,'xiR',OutDir)
call saving3d(XD, 'XD', OutDir)
call saving3d(XN, 'XN', OutDir)

call saving3d(yMat,'yMat',OutDir)

call saving3d(ExpDefProb1,'ExpDefProb1',OutDir)
call saving4d(ExpDefProb,'ExpDefProb',OutDir)

call saving3d(Cond_DefProb1,'Cond_DefProb1',OutDir)
call saving4d(Cond_DefProb,'Cond_DefProb',OutDir)

! grids 
call savingv(aGrid,'aGrid',OutDir)
call savingv(e1Grid,'e1Grid',OutDir)
call savingv(e2Grid,'e2Grid',OutDir)
call savingv(e3Grid,'e3Grid',OutDir)

call savingv(AvgBinSpread,'AvgBinSpread',OutDir)
call savingv(StDevBinSpread,'StDevBinSpread',OutDir)
call savingv(100.0_rk * pDPartition,'pDPartition',OutDir)
call savingv(100.0_rk * pDPartitionPlot,'pDPartitionPlot',OutDir)
call savingv(MassPartition,'MassPartition',OutDir)

call savingm(iTrans,'iTrans',OutDir)

end subroutine SaveSSOutput




















!=====================!
! 					  !
! PRIMITIVE FUNCTIONS !
!					  !
!=====================!
! UTILITY FUNCTION 
elemental function utility(consumption)
real(rk):: consumption, utility
intent(in):: consumption
if (dabs(crra - 1.0_rk).lt.0.000001_rk) then
    utility = dlog(consumption)
else
    utility = consumption**(1.0_rk - crra)
    utility = (utility - 1.0_rk)/(1.0_rk - crra)
end if
end function utility

! COMPUTE LOAN SIZE GIVEN CURRENT AND FUTURE A' INDICES 
pure function LoanSizeInt(ai, api)
integer(ik), intent(in) :: ai, api 
real(rk) :: LoanSizeInt 

if (api .lt. aZeroInd .and. ai .ge. aZeroInd) then
	LoanSizeInt = -aGrid(api)
elseif (api .lt. aZeroInd .and. aGrid(api) - (1.0_rk - phi) * aGrid(ai) .lt. 0.0_rk) then
	LoanSizeInt = (1.0_rk - phi) * aGrid(ai) - aGrid(api)
else 
	LoanSizeInt = 0.0_rk 
end if
end function LoanSizeInt 

! MAP INDEX FROM FULL AGRID ONTO AGRIDPOS
elemental function aiPositive(ai)
integer(ik), intent(in) :: ai
integer(ik) :: aiPositive
aiPositive = ai - naNeg 
end function aiPositive 

! MAP INDEX FROM AGRIDPOS ONTO FULL AGRID
elemental function aiFull(aiPos)
integer(ik), intent(in) :: aiPos
integer(ik) :: aiFull
aiFull = aiPos + naNeg 
end function aiFull 


! DISTRIBUTION FOR THE XI SHOCKS IN DELINQUENCY 
! PDF 
elemental function expPDF(x, mean)
real(rk) :: expPDF, x, mean
intent(in) :: x, mean
expPDF = exp(-x/mean) / mean 
end function expPDF
! CDF 
elemental function expCDF(x,mean)
real(rk) :: expCDF, x, mean
intent(in) :: x, mean
expCDF = 1.0_rk - exp(-x/mean)
end function expCDF 
! inverse CDF 
elemental function expInverseCDF(x,mean)
real(rk) :: expInverseCDF, x, mean
intent(in) :: x, mean
expInverseCDF = -mean * log(1.0_rk - x)
end function expInverseCDF 




!==================================================================!
!												                   !
! COMPUTE CONSUMPTION AND RECOVERY IN CHAPTERS 7 AND 13 BANKRUPTCY !
!												                   !
!==================================================================!
! CHAPTER 7
subroutine Chapter7Details()

integer(ik) :: ai, e1i, e2i, e3i 
real(rk) :: yVal, cVal, aVal, xiVal, xiTildeVal

do e3i = 1,ne3
do e2i = 1,ne2 
do e1i = 1,ne1 
	yVal = yMat(e1i, e2i, e3i)
!$omp parallel do private(aVal, cVal, xiTildeVal, xiVal)
do ai = 1,naNeg
	aVal = aGrid(ai)
	cVal = yVal - psi + xi_7 * aVal
	if (cVal .lt. cMin) then 	! minimum consumption binds
		cVal = cMin
		xiTildeVal = yVal - psi - cMin	! adjust amount of actual recovery
		xiVal = -xiTildeVal / aVal
	else 
		xiTildeVal = -xi_7 * aVal
		xiVal = xi_7
	end if 
	cD7(ai, e1i, e2i, e3i) = cVal
	xiTilde7(ai, e1i, e2i, e3i) = xiTildeVal
	xi7(ai, e1i, e2i, e3i) = xiVal
end do 
!$omp end parallel do 
end do
end do 
end do 

end subroutine Chapter7Details




! CHAPTER 13
subroutine Chapter13Details()

integer(ik) :: e1i, e2i, e3i, i2, i2m, i3, ai
real(rk) :: tSimr, yVal
integer(ik) :: t, tSimCh13
integer(ik) :: j, nSimCh13 = 10000_ik	! number of paths for each individual to simulate. 
real(rk), allocatable, dimension(:) :: yStream
real(rk), allocatable, dimension(:,:) :: e2Rand, e3Rand
real(rk) :: yBarVal, aVal, lambda, xiTildeVal


! SIMULATE EXPECTED DISCOUNTED EARNINGS OVER TIME HORIZON GIVEN THE STATE 
! simulate for a number of periods until discount factor [(1-theta) / (1+r)] ^ i is sufficiently small 
tSimr = log(qTol) / log((1.0_rk - theta_13) / (1.0_rk + rf))
tSimCh13 = floor(tSimr)
print*, 'tSim', tSimCh13

allocate(e2Rand(nSimCh13, tSimCh13), e3Rand(nSimCh13, tSimCh13))
call random_number(e2Rand)
call random_number(e3Rand)

allocate(yStream(nSimCh13))

do e2i = 1,ne2 
do e1i = 1,ne1
	!$omp parallel do private(i2, i2m, i3, yVal)
	do j = 1,nSimCh13
		i2m = e2i	! every individual being simulated starts with the current assumed e2
		yVal = 0.0_rk
		do t = 1,tSimCh13
			i2 = count(e2Rand(j,t) .gt.  cum_e2Trans(i2m,:)) + 1_ik
			i3 = count(e3Rand(j,t) .gt.  cum_e3Prob) + 1_ik
			yVal = yVal + yMat(e1i, i2, i3) * ((1.0_rk - theta_13) / (1.0_rk + rf)) ** real(t)
			i2m = i2
		end do
		yStream(j) = yVal
	end do 
	!$omp end parallel do
	yBarMat(e1i, e2i) = sum(yStream / (1.0_rk - theta_13)) / real(nSimCh13)	! see footnote in Ch 13 section
end do 
end do 


! GIVEN DISCOUNT EXPECTED AVERAGE EARNINGS, COMPUTE THE GARNISHMENT RATE ASSOCIATED WITH EACH STATE. 
do e3i = 1,ne3
do e2i = 1,ne2
do e1i = 1,ne1
	yBarVal = yBarMat(e1i, e1i)
	cD13(e1i, e2i, e3i) = yMat(e1i, e2i, e3i) - psi
	!$omp parallel do private(aVal, lambda, xiTildeVal)
	do ai = 1,naNeg 
		aVal = aGrid(ai)
		lambda = min(-(1.0_rk + rf) * aVal * xi_13 / yBarVal, lambdaBar)
		xiTildeVal = lambda * yBarVal 
		call interp_1d(lpiMat(ai, e1i, e2i, e3i), lpwMat(ai, e1i, e2i, e3i), lambda, nl, lambdaGrid)
		lambda13(ai, e1i, e2i, e3i) = lambda 
		xiTilde13(ai, e1i, e2i, e3i) = xiTildeVal 
		xi13(ai, e1i, e2i, e3i) = xiTildeVal / (-aVal) 
	end do
	!$omp end parallel do 
end do 
end do 
end do 

deallocate(e2Rand, e3Rand, yStream)

end subroutine Chapter13Details







!=================================!
!                                 !
! SIMULATE A PANEL FROM THE MODEL !
!                                 !
!=================================!
! this subroutine draws initial conditions and sets up the panel dimension of the simulation, calling the individual subroutine below 
subroutine PanelSimulation(nSim, tSim, delta_n)
integer(ik), intent(in) :: nSim, tSim 
real(rk), intent(in), dimension(na, ni, ne1) :: delta_n

integer(ik) :: ai, i, e1i, li, j, n, t, nj0, nj7, nj13, aiPos
integer(ik), dimension(nSim) :: aiInit, iInit, e1iInit, liInit
integer(ik), allocatable, dimension(:) :: ai0Vec, i0Vec, e1i0Vec, ai7Vec, i7Vec, e1i7Vec, ai13Vec, i13Vec, e1i13Vec, li13Vec
integer(ik), dimension(nSim) :: flag0 	! initial flag status 
real(rk), allocatable, dimension(:) :: mu0Vec, mu7Vec, mu13Vec, mu0CDF, mu7CDF, mu13CDF
real(rk), dimension(3) :: flagCDF 
real(rk), dimension(nSim) :: InitialRand, flagRand, DefProbSim
real(rk), dimension(nSim, tSim) :: iRand, aPrimeRand, DefRand, DefTypeRand, lambdaRand, FlagRemovalRand
real(rk), dimension(nSim, tSim, 2) :: DischargeRand

real(rk), dimension(nSim, tSim) :: aSim, ySim, DefInd, Ch7Ind, Ch13Ind, RInd, SpreadSim


!-------------------------!
! DRAW INITIAL CONDITIONS !
!-------------------------!
! compute total mass for each flag status, and assign each person in the simulation an initial flag status 
print*, 'flag CDF'
flagCDF(1) = sum(mu0)
flagCDF(2) = flagCDF(1) + sum(mu7)
flagCDF(3) = flagCDF(2) + sum(mu13)
print*, flagCDF
call random_number(flagRand)
do n = 1,nSim
	flag0(n) = count(flagRand(n) .gt.  flagCDF) + 1_ik	! subsequent draws come from transition matrix 
end do

print*, maxval(flag0), minval(flag0)

! set up CDF implied by each stationary distribution 
! -- no flag 
nj0 = na * ni * ne1 
allocate(mu0Vec(nj0), ai0Vec(nj0), i0Vec(nj0), e1i0Vec(nj0))
do e1i = 1,ne1 
do i = 1,ni
do ai = 1,na
	j = (ai - 1_ik) * ni * ne1 + (i - 1_ik) * ne1 + e1i 
	mu0Vec(j) = mu0(ai, i, e1i)
	ai0Vec(j) = ai
	i0Vec(j) = i
	e1i0Vec(j) = e1i 
end do 
end do
end do 
mu0Vec = mu0Vec / sum(mu0Vec)

allocate(mu0CDF(nj0))
call CreateCDF(mu0CDF, mu0Vec, nj0)
deallocate(mu0Vec)


! -- Chapter 7
nj7 = naPos * ni * ne1 
allocate(mu7Vec(nj7), ai7Vec(nj7), i7Vec(nj7), e1i7Vec(nj7))
do e1i = 1,ne1 
do i = 1,ni
do ai = 1,naPos
	j = (ai - 1_ik) * ni * ne1 + (i - 1_ik) * ne1 + e1i 
	mu7Vec(j) = mu7(ai, i, e1i)
	ai7Vec(j) = ai
	i7Vec(j) = i
	e1i7Vec(j) = e1i 
end do 
end do
end do 
mu7Vec = mu7Vec / sum(mu7Vec)

allocate(mu7CDF(nj7))
call CreateCDF(mu7CDF, mu7Vec, nj7)
deallocate(mu7Vec)


! -- Chapter 13 
nj13 = naPos * ni * ne1 * nl
allocate(mu13Vec(nj13), ai13Vec(nj13), i13Vec(nj13), e1i13Vec(nj13), li13Vec(nj13))
do li = 1,nl
do e1i = 1,ne1 
do i = 1,ni
do ai = 1,naPos
	j = (ai - 1_ik) * ni * ne1 * nl + (i - 1_ik) * ne1 * nl + (e1i - 1_ik) * nl + li 
	mu13Vec(j) = mu13(ai, i, e1i, li)
	ai13Vec(j) = ai
	i13Vec(j) = i
	e1i13Vec(j) = e1i 
	li13Vec(j) = li 
end do 
end do
end do 
end do 
mu13Vec = mu13Vec / sum(mu13Vec)

allocate(mu13CDF(nj13))
call CreateCDF(mu13CDF, mu13Vec, nj13)
deallocate(mu13Vec)


! given flag statuses determined above, draw other parts of initial state using the CDFs just computed 
call random_number(InitialRand)
777 format(10i10)
do n = 1,nSim 
	if (flag0(n) .eq. 1_ik) then 			! no flag 
		j = count(InitialRand(n) .gt.  mu0CDF) + 1_ik
		aiInit(n) = ai0Vec(j)
		iInit(n) = i0Vec(j)
		e1iInit(n) = e1i0Vec(j)
		liInit(n) = 1_ik					! pure placeholder 
	elseif (flag0(n) .eq. 2_ik) then 		! chapter 7 flag 
		j = count(InitialRand(n) .gt.  mu7CDF) + 1_ik
		aiInit(n) = aiFull(ai7Vec(j))		! keep everything on the full scale to keep things straight 
		iInit(n) = i7Vec(j)
		e1iInit(n) = e1i7Vec(j)
		liInit(n) = 1_ik					! pure placeholder 
	else 									! chapter 13 flag 
		j = count(InitialRand(n) .gt.  mu13CDF) + 1_ik
		aiInit(n) = aiFull(ai13Vec(j))		! keep everything on the full scale to keep things straight 
		iInit(n) = i13Vec(j)
		e1iInit(n) = e1i13Vec(j)
		liInit(n) = li13Vec(j)
	end if 
end do 

print*, ' '



!-----------------------------------------------------------------------------!
! CONSTRUCT CUMULATIVE VERSIONS OF ALL DECISION RULES AND TRANSITION MATRICES !
!-----------------------------------------------------------------------------!
allocate(cum_iTrans(ni, ni))
! exogenous idiosyncratic states 
do i = 1,ni 
	call CreateCDF(cum_iTrans(i,:), iTrans(i,:), ni)
end do

! discharge shocks 
allocate(cum_xiProb(nXi))
call CreateCDF(cum_xiProb, xiProb, nXi)

! decision rules 
allocate(cum_pA0(na, na, ni, ne1))
allocate(cum_pA7(naPos, naPos, ni, ne1))
allocate(cum_pA13(naPos, naPos, ni, ne1, nl))
do e1i = 1,ne1
do i = 1,ni
do ai = 1,na 
	call CreateCDF(cum_pA0(:, ai, i, e1i), pA0(:, ai, i, e1i), na)
	if (ai .ge. aZeroInd) then 			! determine which default type conditional on default
		aiPos = aiPositive(ai)
		call CreateCDF(cum_pA7(:, aiPos, i, e1i), pA7(:, aiPos, i, e1i), naPos)
		do li = 1,nl
			call CreateCDF(cum_pA13(:, aiPos, i, e1i, li), pA13(:, aiPos, i, e1i, li), naPos)
		end do 
	end if 
end do
end do 
end do 


!---------------------------------------------!
! LOOP OVER INDIVIDUALS TO SIMULATE THE PANEL !
!---------------------------------------------!
! DRAW RANDOM NUMBERS
call random_number(DefRand)
call random_number(DefTypeRand)
call random_number(lambdaRand)
call random_number(aPrimeRand)
call random_number(iRand)
call random_number(DischargeRand)
call random_number(FlagRemovalRand)

! LOOP OVER INDIVIDUALS 
print*, 'simulations'
!$omp parallel do
do n = 1,nSim 
	DefProbSim(n) = delta_n(aiInit(n), iInit(n), e1iInit(n))
	call SimulateIndividual(DefInd(n,:), Ch7Ind(n,:), Ch13Ind(n,:), RInd(n,:), SpreadSim(n,:), aSim(n,:), ySim(n,:), &
		& tSim, flag0(n), aiInit(n), iInit(n), e1iInit(n), liInit(n), &
		& DefRand(n,:), DefTypeRand(n,:), lambdaRand(n,:), aPrimeRand(n,:), iRand(n,:), FlagRemovalRand(n,:), DischargeRand(n,:,:))
end do 
!$omp end parallel do 


!-----------------!
! COMPUTE MOMENTS !
!-----------------!
900 format(10a10)
print 900, 't','BK rate','Ch 7 %', 'Ch 13 %','avg a','avg y'
901 format(i10, 9f10.4)
do t = 1,tSim
	print 901, t, &
		& sum(Ch7Ind(:,t) + Ch13Ind(:,t)) / real(nSim), &
		& sum(Ch7Ind(:,t)) / sum(Ch7Ind(:,t) + Ch13Ind(:,t)), &
		& sum(Ch13Ind(:,t)) / sum(Ch7Ind(:,t) + Ch13Ind(:,t)), &
		& sum(aSim(:,t)) / real(nSim), &
		& sum(ySim(:,t)) / real(nSim)
end do 

print*, 'saving panel output'
call savings(dble(tSim),'tSim',PanelDir)
call savings(dble(nSim),'nSim',PanelDir)

call savingv(DefProbSim,'DefProbSim',PanelDir)

call savingm(DefInd,'DefInd',PanelDir)
call savingm(Ch7Ind,'Ch7Ind',PanelDir)
call savingm(Ch13Ind,'Ch13Ind',PanelDir)
call savingm(RInd,'RInd',PanelDir)
call savingm(SpreadSim,'SpreadSim',PanelDir)
call savingm(aSim,'aSim',PanelDir)
call savingm(ySim,'ySim',PanelDir)

end subroutine PanelSimulation 










! this subroutine simulates a given individual for T periods, given all the relevant draws of randomness
subroutine SimulateIndividual(DefInd, Ch7Ind, Ch13Ind, RInd, SpreadSim, aSim, ySim, &
	& tSim, flag0, ai0, i0, e1i, li0, DefRand, DefTypeRand, lambdaRand, aPrimeRand, iRand, FlagRemovalRand, DischargeRand)
	
integer(ik), intent(in) :: tSim, flag0, e1i, ai0, i0, li0
real(rk), intent(in), dimension(tSim) :: DefRand, DefTypeRand, lambdaRand, aPrimeRand, iRand, FlagRemovalRand
real(rk), intent(in), dimension(tSim, 2) :: DischargeRand

real(rk), intent(out), dimension(tSim) :: DefInd, Ch7Ind, Ch13Ind, RInd, SpreadSim, aSim, ySim

integer(ik) :: t, j, flag, flagp, i, ip, ai, api, e2i, e3i, aiPos, apiPos, li, lpi
real(rk) :: apw, lpw 

! INITIALIZATIONS 
DefInd = 0.0_rk 
Ch7Ind = 0.0_rk
Ch13Ind = 0.0_rk 
RInd = 0.0_rk 
SpreadSim = 0.0_rk 

! SET UP FIRST ITERATION
flag = flag0  
ai = ai0
i = i0
li = li0

! LOOP THROUGH PERIODS
do t = 1,tSim 
	e2i = e2Map(i)
	e3i = e3Map(i)
	
	aSim(t) = aGrid(ai)
	ySim(t) = yMat(e1i, e2i, e3i)
	
	! NO FLAG 
	if (flag .eq. 1_ik) then
		! If in debt, have to consider default options 
		if (ai .lt. aZeroInd) then 
		
			!---------!
			! DEFAULT !
			!---------!
			if (DefRand(t) .le. pD(ai, i, e1i)) then
				DefInd(t) = 1.0_rk 
				if (DefTypeRand(t) .le. pD7(ai, i, e1i)) then 							! CH 7
					Ch7Ind(t) = 1.0_rk 
					api = aZeroInd 
					flagp = 2_ik
					lpi = 1_ik 
				elseif (DefTypeRand(t) .le. pD7(ai, i, e1i) + pD13(ai, i, e1i)) then 	! CH 13
					Ch13Ind(t) = 1.0_rk 
					api = aZeroInd 
					flagp = 3_ik
					if (lambdaRand(t) .le. lpwMat(ai, e1i, e2i, e3i)) then 
						lpi = lpiMat(ai, e1i, e2i, e3i)
					else
						lpi = lpiMat(ai, e1i, e2i, e3i) + 1_ik
					end if 
				else 																	! Reorganization
					RInd(t) = 1.0_rk 
					flagp = 1_ik
					lpi = 1_ik
					j = count(DischargeRand(t,1) .gt. cum_xiProb) + 1_ik
					call interp_1d(api, apw, (1.0_rk - xiGrid(j)) * aGrid(ai), na, aGrid)
					if (DischargeRand(t,2) .le. apw) then 
						api = api
					else 
						api = api + 1_ik
					end if 
				end if 
				
			!------------!
			! NO DEFAULT !
			!------------!
			else 
				api = count(aPrimeRand(t) .gt. cum_pA0(:,ai,i,e1i)) + 1_ik 
				if (api .lt. aZeroInd) then 
					SpreadSim(t) = spreads(api, ai, i, e1i)
				end if 
				lpi = 1_ik
				flagp = 1_ik
			end if 
		
		! If not in debt, only consider savings / borrowing 
		else
			api = count(aPrimeRand(t) .gt. cum_pA0(:,ai,i,e1i)) + 1_ik 
			if (api .lt. aZeroInd) then 
				SpreadSim(t) = spreads(api, ai, i, e1i)
			end if 
			lpi = 1_ik
			flagp = 1_ik
		end if 
		
	! CH 7 FLAG -- only consider saving 
	elseif (flag .eq. 2_ik) then
		aiPos = aiPositive(ai)
		apiPos = count(aPrimeRand(t) .gt. cum_pA7(:,aiPos,i,e1i)) + 1_ik 
		api = aiFull(apiPos)
		lpi = 1_ik 
		if (FlagRemovalRand(t) .le. theta_7) then 
			flagp = 1_ik
		else 
			flagp = 2_ik 
		end if 
	! CH 13 FLAG -- only consider saving 
	else 
		aiPos = aiPositive(ai)
		apiPos = count(aPrimeRand(t) .gt. cum_pA13(:,aiPos,i,e1i,li)) + 1_ik 
		api = aiFull(apiPos)
		lpi = li
		if (FlagRemovalRand(t) .lt. theta_13) then 
			flagp = 1_ik
		else 
			flagp = 3_ik 
		end if 
	end if 
	
	
	!---------------------------------------------!
	! SET UP THE NEXT ITERATION (if not the last) !
	!---------------------------------------------!
	if (t .lt. tSim) then 
		! draw exogenous idiosyncratic shocks
		ip = count(iRand(t+1) .gt.  cum_iTrans(i,:)) + 1_ik
		
		! update the state 
		flag = flagp
		li = lpi 
		ai = api 
		i = ip
	end if 
end do 

end subroutine SimulateIndividual





subroutine CreateCDF(CDF, PDF, n)
integer(ik), intent(in) :: n
real(rk), intent(in), dimension(n) :: PDF 
real(rk), intent(out), dimension(n) :: CDF 
integer(ik) :: i 
CDF(1) = PDF(1)
do i = 2,n
	CDF(i) = CDF(i-1) + PDF(i)
end do 
CDF = CDF / sum(PDF)
end subroutine CreateCDF



end module SS_JPEM_Module